SM472, Computational Fourier Transforms u 

Class notes for a Capstone course taught Spring 2006-2007. The text 
for the course is the book Fast Fourier transform by James Walker |Wlj . 
Examples using SAGE [S] illustrate the computations are given throughout. 
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Application of FS, FTs, and DFTs (including DCTs and FFTs), are found 

in 

• the heat equation and the wave equation in physics, 

• image compression, 

• analog-to-digital conversion, 

• spectroscopy (in medicine, for example), 

• passive sonar (via frequency analysis), 

• multiphcation of very large integers (useful in cryptography), 

• statistics, 

among many others. 

In fact. Patent number 5835392 is "Method for performing complex fast 
fourier transforms (FFT's)", and Patent number 6609140 is "Methods and 
apparatus for fast fourier transforms." These are just two of the many patents 
which refer to FFTs. Patent number 4025769 is "Apparatus for convoluting 
complex functions expressed as fourier series" and Patent number 5596661 
is "Monolithic optical waveguide filters based on Fourier expansion." These 
are just two of the many patents which refer to FSs. You can look these up 
in Google's Patent search engine. 
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http : //www . google . com/ptshp?hl=en&tab=wt&q=^ 



to read more details. You'll see related patents on the same page. Try 
"discrete cosine transform" or "digital filters" to see other devices which are 
related to the topic of this course. 



1 Background on Fourier transforms 

This section is not, strictly speaking, needed for the introduction of discrete 
Fourier transforms, it helps put things in the right conceptual framework. 



1.1 Motivation 

First, let's start with some motivation for the form of some of the defini- 
tions. Let's start with something you are very familiar with: power series 
^^Q/(t)a;*, where we are writing the coefficients f{t), t G N, of the power 
series in a functional notation instead of using subscripts. Usually, if /(O), 
/(I), /(2), ... is a given sequence then the power series X^t^o /(^)^* called 
the generating function of the sequenc^ 

The continuous analog of a sum is an integral: 

/ fit)x'dt. 

t=o -^0 
Replacing x by e"**, the map 

/ f{t)x'dt= / f{t)e-''dt, 



is the Laplace transform. Replacing xhy e the map 

f^oo poo 

/ f{t)x'dt= / f{t)e~'^'ydt, 



is essentially the Fourier transform (the definition below includes a factor 
for convenience). In other words, the Laplace transform and the Fourier 
transform are both continuous analogs of power series. 



■^In fact, since power series are so well understood, often times one studies the generating 
function to better understand a given sequence. 
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To go a little further with this analogy, suppose we have two power series 
A{x) = Yl'jLo^U)^'' B{x) = '^'kLoB{k)x''. The product is given by 

A{x)B{x) = (Er=o«o>^)(Er=o^(^)^') 

= (a(0) + a{l)x + a{2)x^ + ...){b{0) + b{l)x + b{2)x^ + ...) 

= a(0)6(0) + (a(0)6(l) + a(l)6(0))x + (a(2)6(0) + a(l)6(l) + a{0)b{2))x^ + ... 

where 

m 

Cm= ^(•^')^(^) = - j)- 

j+k=m j=0 

This last expression is referred to as the Dirichle^ convolution of the a/s 
and 6fc's. Likewise, if F{s) = f{t)e~^^dt and G{s) = g{t)e~'^^ dt then 

F{s)G{s)= / {f * g){t)e~^Ut, 
Jo 

where 

(/*^)(t)= f f{z)g{t-z)dz. 
Jo 

The above integral, f{t)e~'^^^y dt, is over < t < oo. Equivalently, it 
is an integral over M, but restricted to functions which vanish off (0, oo). If 
you replace these functions supported on (0, oo) by any function, then we are 
lead to the transform 

/oo 
f{t)e-'''ydt. 
'OO 

This is discussed in the next section. 

1.2 The Fourier transform on M 

If / is any function on the real line M whose absolute value is integrable 

||/||i= / \f{x)\dx <oo 
Jr 

■^Pronounced "Dear-ish-lay" . 
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then we say / is L^, written / £ L^(]R). In this case, we define the Fourier 
transform of / by 




This is a bounded function, written / e L°°(]R), 

I I/I loo = sup 1/(2/) I < 00. 

It can be shown that F (extends to a well-defined operator which) sends 
any square-integrable function to another square-integrable function. In 
other words, 

F : L2(M) ^ 

where L^(M) denotes the vector space (actually, a Hilbert space) of functions 
for which the L^-norm is finite: 

ll/l|2=(/ \f{x)\'dxy/'<oo. 

Jr 

Example 1 Here is an example. Let fo{x) = e~'^^ . We have 

^ JIm(z)=y ^ 

where z = x+iy the last equality follows from the Residue Theorem from complex analysis. 
Now, this is 

= e-y\ 

This last line depends on 




(Here's the quick proof of this: 
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Now take square-roots.) 

This proves that this example satisfies 



F{fo){y) - /o(y). 



(1) 



In other words, that the operator F has eigenvector /o and eigenvalue A = 1. 
Exercises: 

1. Verify that F{f'){y) = 2iyF{f){y), for all / G ^^(R) n C\R). 
(Hint: Integrate by parts.) 

2. If fi{x) = xe~^^ then /i is an eigenvector of the Fourier transform with 
eigenvalue A = —i, i.e., = —ifi- 

(Hint: Use ((H) and the previous exercise.) 



3. Let 




and show that f{y) 



sinj-Ky) 
Try 



(= sinc(?/)). 



4. For a < 0, let 




and show that f{y) 



{2TT){(iy-a) 



. Similarly, let 




and show that f{y) 



(27r)((iy+a) " 



-1 



5. For a < 0, let /(x) 



g27ra|x| show that /(?/) 



W{(s/2+a2)- 



a 
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The following result spells out the main properties of the FT 
Theorem 2 Let f,g e L\R), with f G L\R) and xf{x) e 
• Linearity: 



• Scaling: 



and 



Shifting: 



Modulation: 



Inversion: 



Differentiation: 



and 



af + bg i — > af + hg, 



f{x/a) I — > af{ay), 



f{ax) ^ a ^f{y/a), 



fix-a)^e-''^^^yf{y), 



e2--/(a;)^/(y-a), 



/oo 
m 
■oo 



e'^^*^ dt. 



fix) ^ 2myf{y), 



Parseval: For f,g e L^(M) n ^^(M), we have 

f{x)g{x)dx= / f{y)g{y)dy, 
Jm. 



and in particular, 



\fix)\'dx= / \f{y)\'dy 



This can be found in Chapter 5 of Walker |Wlj . 



1.3 Application to Laplace's PDE 

We shall use the FT to find a function w{x,t) satisfying 



9x2 -r Q^2 



0, 



w{x,0) = f{x). 

This represents the steady state temperature of a semi-infinite plate, in the 
shape of the upper half-plane, where the boundary (the x-axis) is kept heated 
to the temperature /(x) at the point (x, 0) for all time. 

solution: Assume that for each y, w{x,y) is in L^(M), as a function of x. 

Let 

w{u,y)= I w(a;, ?/)e"^™"(ia;, 

so 

w(m,0)= f w;(x,0)e-2"^^^da;,= f /(x)e-2™" dx = /(m), 
by the initial condition, and 



w{x,y)= / w{u,y)e '^du 



by the inversion formula. 
By (j2I), we have 



w^^{x,y) 



-(27rM)^u;(M,y)e'™"rfM 



and 



Wyy{x,y)= I Wyy{u,y)e^''"''du, 



(2) 
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so 

= Wxx{x, y) + Wyy{x, y) = {-{2nufw{u, y) + Wyy{u, 2/))e^''*^" du. 

Jr 

Therefore, by Laplace's equation, 

-{2'Ku)^w{u, y) + Wyyiu, y) = 0, 

or 

k"iy) - {27:ufk{y) = 0, 
where k{y) = w{u,y). This means 

w{u,y) = k{y) = + C2e2'^l"l^ 



for some constants ci,C2. Let C2 = (because it works, that's why!) and 
solve for ci using the initial condition w{u, 0) = f{u), to get 

wiu,y) = fiu)e-'''\-\y. 

Using (j21) again, we have 

w{x,y) = /jg/(u)e-2^l"M^^^"rfM 

= IrUr e-^^l^l^'e^'^*^^-^)" du)f{z) dz 

which is the convolution of /(x) with Ky{x) = ^ y2^^2 - This function w{x, y) = 
(/ * Ky){x) is the desired solution. 

1.4 Application to Schodinger's PDE 

Assume / G L\R)nL^{R). 

We shall use the FT to find a function ilj{x,t) satisfying 

dx^ dt ' 

^/'(x,0)=/(a;). 
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where (for convenience) o,^ = and h is Planck's constant. The function 
'ilj{x,t)\'^, where ip is the solution (or "wave function") of the above PDE, 
represents (under certain conditions on /) the probability density function 
of a free particle on a semi-infinite plate, in the shape of the upper half-plane, 
where the distribution on boundary (the x-axis) is determined by f{x). 
solution: Assume that for each t, ip{x^t) is in //-"^(M), as a function of x. 

Let 

i){u,t) = I ^(x,t)e-2"""rfa;, 

SO 

^(m,0)= I V^(x,0)e"2"""rfx,= f /(x)e"2"""rfx = /(m), 
Jr Jr 

by the initial condition, and 

i){x,t) = [ V^(M,t)e2™"du, (3) 
Jr 

by the inversion formula. By Q, we have 

^..{x,t) = /-(27rM)2^(M,t)e2--"rfM, 

and 

Jr 

so 

= a^^,,{x, t) - ^<(x, t)= [ {-a\2nuftijiu, t) - ^^(m, t))e2'^*"" du. 

Jr 

By Schodinger's PDE, 

-a^(27rM)2^(M, t) - ^tiu, t) = 0, 

or 

k'{t) = -a'{2T:ufk{t), 



10 



where k{t) = ■?/'(«, t). This means 



^^J{u,t) = k{t) = coe""' 

for some constant Cq. Using the above initial condition to sove for cq, we 
obtain ij{u,t) = f{u)e~°' *. Using again, we have 

= ImHtix - z)f{z) dz, 

which is the convolution of f{x) with Ht{x) = Jj^ e~"^*^^™^^*e^'^*'*" du. This 
function w{x,y) = (/ * Ht){x) is the desired solution. 
By the Plancheral formula, 

lR\i^ix,t)? dx = f^\ij{u,t)\^du 

= /rI/(w)Nw 
= 1^11(^)1^ dx 
= /Rl^(w,0)|2da;, 

by the initial condition. In particular, if / has L^-norm equal to 1 then so 
does ■ilj{x,t), for each t>0. 

2 Fourier series 

History: Fourier series were discovered by J. Fourier, a Frenchman who 
was a mathematician among other things. In fact, Fourier was Napolean's 
scientific advisor during France's invasion of Egypt in the late 1800's. When 
Napolean returned to France, he "elected" (i.e., appointed) Fourier to be a 
Prefect - basically an important administrative post where he oversaw some 
large construction projects, such as highway constructions. It was during this 
time when Fourier worked on the theory of heat on the side. His solution to 
the heat equation is basically what we teach in the last few weeks of your 
differential equations course. The exception being that our understanding of 
Fourier series now is much better than what was known in the early 1800's. 
For example, Fourier did not know of the integral formulas (jlj), ^ for the 
Fourier coefficients given below. 
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Motivation: Fourier series, since series, and cosine series are all expan- 
sions for a function f{x) in terms of trigonometric functions, much in the 
same way that a Taylor series is an expansion in terms of power functions. 
Both Fourier and Taylor series can be used to approximate a sufficient ;y 
"nice" function f{x). There are at least three important differences between 
the two types of series. (1) For a function to have a Taylor series it must 
be different iabl^ whereas for a Fourier series it does not even have to be 
continuous. (2) Another difference is that the Taylor series is typically not 
periodic (though it can be in some cases), whereas a Fourier series is always 
periodic. (3) Finally, the Taylor series (when it converges) always converges 
to the function f{x), but the Fourier series may not (see Diri chiefs theorem 
below for a more precise description of what happens). 

Given a "nice" periodic function f{x) of period P, there are with n > 
and bn with n > 1 such that f{x) has ("real") Fourier series 

mr, //.N/ X ^0 v-^r ,2'Knx. , . .2Tinx.^ 
FS^[f)[x) = Y + 2^[a„cos(— ^)+6„sm(— ^)]. 

n=l 

These Fourier series coefficients are given by 

2 /"^ ,27mx 

an = 

and 



^^=pI /(3;)cos(— ^)rfx, (4) 



/(x) sin( — ——)dx. (5) 



" p 

What does "nice" mean? When does the series converge? When it does 
converge, what is the relationship between / and FS'k(/)? 

Definition 3 A function f{x) on a finite interval [a,b] is said to be of bounded 
variation if there is a constanct C > such that for any partition xq = a < xi < 
... < Xn-i < Xn = b of the interval (a, 6), we have 

n-l 



^|/(x,+i)-/(x,)| <C. 



4 



i=0 

Remember the formula for the n-th Taylor series coefficient centered at a; = a 

n! 



_ /'"'(a) v 
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Another characterization states that the functions of bounded variation on a 
closed interval are exactly those functions which can be written as a difference 
g — h, where both g and h are monotone. 

Let BV{[0, P]) denote the vector space of functions on [0, P] of bounded vari- 
ation. 

Since any monotone function is Riemann integrable, so is any function 
of bounded variation. The map FS is therefore well-defined on elements of 
BVi[0,P]). 

Likewise, given a "nice" periodic function f{x) of period P, there are a„ 
with n > and 6„ with n> 1 such that f{x) has ("complex") Fourier series 



FS{f){x)= J2 Cne— . 



These Fourier series coefficients are given by 



cn = -^^ /(x)e-^™^/^rfa;, (6) 
for n G Z. 

Theorem 4 /// G BV{[0, P]) and if f has a Fourier series representation 



oo 

2-Kinx 



n=— oo 

which converges absolutely then ^ holds. 
proof: By hypothesis, 

„p „p oo 

•^0 ^0 rn=-oo 

Interchanging the sum and the integral, this is 



oo „p 

V c,„ / e^-^^'"-")^/^ dx = Pc„ 



m=— oo 



□ 

In the above proof, we have used the fact that {e^'^*"'^/^}„gz forms a 
sequence of periodic orthogonal functions on the interval (0, P). 
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Example 5 Given any complex Fourier series f{t) = x{t) + iy{t), you can plot the 
parametric curve {x{t),y(t), for < t < P. A. Robert |Q classifies the complex 
Fourier series whose parametric plot is a regular polygon. The n-gons are all 
obtained by transforming in various ways (such as translations and reflections) the 
basic Fourier series 

f{t) = ^(1 + A:n)-2g^i(i+fcn)t. 

Here are some examples. In the first row of the table of plots below, we plot 
the partial sum 

(1 + A:n)-2e'^*(i+^")* 

\k\<N 

with N = 3, n = 3, then = 10, n = 3. The next line is with = 3, n = 4, then 
A^ = 10, n = 4; the last line is with A^ = 3, n = 5, then A^ = 10, n = 5. 
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2.S 5 75 ■ 10 







r-t%it "f«4 t 





Here is the SAGE code which produced these examples: 

sage : N = 3 

sage: L = [l+5*k for k in range(-N,N)] 

sage: f = lambda t: suin( [10*n" (-2)*exp(n*pi*I*t) for n in L] ) 

sage: pts = [[real(f (i/100)) ,imag(f (i/100))] for i in range(200)] 

sage: show(list_plot(pts)) 
sage : 
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sage: N = 10 

sage: L = [l+5*k for k in range (-N,N)] 

sage: f = lambda t: suinC [10*n" (-2)*exp(n*pi*I*t) for n in L] ) 

sage: pts = [[real(f (i/100)) ,imag(f (i/100))] for i in range(200)] 

sage: show(list_plot(pts)) 

This code produces the pentagon pictures above. Changing the terms in the 
partial FS obviously changes the graph. For instance, replacing L = [l+5*k for 
k in range (-N,N)] by L = [2+5*k+k**2 for k in range (-N,N)] produces a 
figure looking like a goldfish! 



The relationship between the real and the complex Fourier series is as 
follows: In FS{f), replace the cosine and sine terms by 



and 



. 27Tnx , 
cos(— ^) 

. , 2nnx , 



^■ninxjP ^—2mnxjP 



sm 



P 



2i 



This implies FSM,{f ) = FS{f), where 



Cn = < 



ap 
2 ' 



2 2i ' 



n = 0, 
n > 0, 



. 2 2i ' "' ^ ^■ 

In particular, FSM.{f){x) converges if and only if FS{f){x) does. 

Let C'^{P) denote the vector space of all functions which are m-times 
continuously differentiable and have period P. 

Lemma 6 // / G C"^{P) then there is a constant A = A{f) > such that 
\cn\ < ^/l^^'l ™; for all n 7^ 0. 



proof: Integrate-by-parts: 
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since /(O) = f{P). 

Now use the "trivial" estimate 

c/(a;)e-2™^/^da;| < P- max \g{x)\ 

0<x<P 

to get 

maxo<a;<p [/'(x)! 1 

The proof is completed using mathematical induction (or simply repeat this 
process m — 1 more times). □ 

Let £°°(Z) denote the space of vectors {xn)n£Z whose coordinates are 
bounded and let 

||(a;„)n6z|| = ||(a;„.)„ez||oo = max \xn\. 
In particular, it follows that the Fourier transform defines a linear mapping 

FS : C\P) ^r{%), 

given by sending a continuous periodic function / G C*^(P) to its sequence 
of Fourier series coefficients (c„(/))„gz ^ ^°^(Z). 

Theorem 7 (Jordan) If f e BV{[0, P]) then FS{f){x) converges to /(^-+)+/(^-) , 
proof: See §13.232 in [!].□ 

Perhaps more familiar is the following result, which was covered in your 
Differential Equations course. 

Theorem 8 (Dirichlet) If f has a finite number of maxima, a finite number 
of minima, and a finite number of discontinuities on the interval (0,P) then 
FS{f){x) converges to iIi±l±Z(2-ii . 

proof: See §13.232 in [TJ.D 

Example 9 Let P = 2 and let f{x) = x for < x < 2. This is sometimes 
referred to as the "sawtooth function", due to the shape of its graph. The 
Fourier coefficients are given by 
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1 ^g-Trmx^^^ lr/^_Wg-7rmx 







when n 0, and 



1 / W ^ 



[—TiinY 



x=2 



x=Q —Trm 



x=2 



x=0 



1. 



The Fourier series 



nez 



neZ, n^^O 



c?oes not converge absolutely. 



Example 10 Let P ^ 1 and let f{x) = 10x2(1 - xY for < x < 1. This 
function belongs to C^(l). 




Figure 1: Graph of f{x), < x < 1. 
The Fourier coefficients are given by 



' lox^i-xfe-^-^ dx = io( ^ ^ - tiir ^^h , 



4 i n^TT^ 



50 Cn = when n ^ 0. When n — 0, we have 



4 i n^n^ 



Co 



= [ 10x'^{l-xfdx^l/3. 
Jo 
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In this case, the Fourier series 



does converge absolutely. In fact, it can even be differentiated term-by-term. 
The SAGE commands 

sage: f = maximaC"10*x**2* Cl-x) **2*expC-2*Pi*i*x*n) " ) 
sage: f . integral C 'x' , 0, 1) 

10*CCi"2*n"2*Pi"2 - 3*i*n*Pi + 3)/ C4*i"5*n"5*Pi"5) - Ci"2*n"2*Pi"2 + 3*i*n*Pi + 3) *°/.e"- C2*i*n*Pi) /C4*i"5*n"5*Pi"5) ) 



were used to compute the Fourier coefficients above. 



2.1 Sine series and cosine series 

Recall, to have a Fourier series you must be given two things: (1) a period 
P = 2L, (2) a function f{x) defined on an interval of length 2L, usually we 
take —L < x < L (but sometimes < x < 2L is used instead). The Fourier 
series of /(x) with period 2L is 

oo 

«o v^r ,mrx. , . .mix.. 

/ W ~ Y + 2^[«n C0S(-^) + hn Sm(-^)J, 
n=l 

where and 6„ are as in (jH), 

First, we discuss cosine series. To have a cosine series you must be given 
two things: (1) a period P = 2L, (2) a function f{x) defined on the interval 
of length L, < X < L. The cosine series of f{x) with period 2L is 



where a„ is given by 



fix) ~ y + 2^ «n cos(-^), 

n=l 



2 ,mix 



L 



cos{——)f{x) dx. 







The cosine series of f{x) is exactly the same as the Fourier series of the even 
extension of f{x), defined by 



feven (2^) 



f{x), < X < L, 
f{-x), -L < X <0. 
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Next, we define sine series. To have a sine series you must be given two 
things: (1) a "period" P = 2L, (2) a function f{x) defined on the interval of 
length L, < X < L. The sine series of f{x) with period 2L is 



fi^) ~ sin 



nrcx , 

n=\ 



where bn is given by 



bn = [ sin(^^)/(x) dx. 



The sine series of f{x) is exactly the same as the Fourier series of the odd 
extension of /(x), defined by 



fodd{x) 



f{x), < X < L, 
-f{-x), -L < X <0. 

One last definition: the symbol ~ is used above instead of = because of 
the fact that was pointed out above: the Fourier series may not converge to 
f{x) at every point (recall Dirichlet's Theorem |SI). 

Example 11 If f{x) = 2 + x, —2 < x < 2, is extended periodically to M with 
period 4 then L = 2. Without even computing the Fourier series, we can evaluate 
the FS using Dirichlet's theorem (Theorem |S1 above). 

Question: Using periodicity and Dirichlet's theorem, find the value that the 
Fourier series of /(x) converges to at x = 1,2,3. (Ans: f{x) is continuous at 1, so 
the FS at X = 1 converges to /(I) — 3 by Dirichlet's theorem. f{x) is not defined at 2. 
It's FS is periodic with period 4, so at a; = 2 the FS converges to ^ o±4 ^ 2. 

f{x) is not defined at 3. It's FS is periodic with period 4, so at a; = 3 the FS converges to 
_ 1+1 _ n 

2 2 ' 

The formulas for and bn enable us to compute the Fourier series coefficients 
ao, an and 6„. These formulas give that the Fourier series of /(x) is 

,./ N , nvr COS (nvr) . .nvrx, 
/ X -4 + 2^-4 sm^ . 

n=l 

The Fourier series approximations to /(x) are 



5o = 2, 5, = 2 + isin(^), 5^ = 2 + 4 ^^^iil^ - 2 ^^^ii^ , ... 

TT 2 TT TT 
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The graphs of each of these functions get closer and closer to the graph of f{x) 
on the interval —2<x<2. For instance, the graph of f{x) and of Ss are given 
below: 




Notice that f{x) is only defined from — 2 < x < 2 yet the Fourier series is not only 
defined everywhere but is periodic with period P = 2L = 4. Also, notice that Ss 
is not a bad approximation to f{x), especially away from its jump discontinuities. 

Example 12 This time, let's consider an example of a cosine series. In this case, 
we take the piecewise constant function f(x) defined on < x < 3 by 



fix) 



1, < X < 2, 

-1, 2 < X < 3. 

We see therefore L = 3. The formula above for the cosine series coefficients gives 
that 
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n=l 

The first few partial sums are 

\/3cos (4 TTx) 

TT 

5*3 = 13 + 2- 

TT TT 

Also, notice that the cosine series approximation Siq is an even function but f{x) 
is not (it's only defined from < x < 3). As before, the more terms in the cosine 
series we take, the better the approximation is, for < x < 3. For instance, the 
graph of f{x) and of are given below: 
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0.5- 



-0.5 




Figure 3: Graph of f{x) and a cosine series approximation of f{x). 

Roughly speaking, the more (everywhere) difFerentiable the function is, the 
faster the Fourier series converges and, therefore, the better the partial sums of 
the Fourier series will approximate f{x). 

Example 13 Finally, let's consider an example of a sine series. In this case, we 
take the piecewise constant function f{x) defined on < x < 3 by the same 
expression we used in the cosine series example above. 

Question: Using periodicity and Dirichlet's theorem, find the value that the 
sine series of f{x) converges to at x = 1, 2, 3. (Ans: f{x) is continuous at 1, so 
the FS at X = 1 converges to /(I) = 1. f{x) is not continuous at 2, so at x = 2 the SS 
converges to = /(-2+)+/(2-) ^ ^ ^ j-^^^ defined at 3. It's SS is 

periodic with period 6, so at a; = 3 the SS converges to ■^°'i'i(^~)+-^°<i<i(^+) = ^^ti = g.) 

The formula above for the sine series coefficients give that 



cos (nir) — 2 cos (I mr) + 1 ,mrx, 

sin — - . 

mr 3 
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The partial sums are 



^2 



^ sin(l/37ra;) ^ ^ sin(|7rx) 



TT 



TT 



S3 = 2 



sin (I TTx) 



TT 



+ 3- 



sm (I irx) 



-4/3 



sin (tt x) 



These partial sums Sfi^ as ?i 



TT TT 

00, converge to their limit about as fast as those 



in the previous example. Instead of taking only 10 terms, this time we take 40. 
Observe from the graph below that the value of the sine series at a:; = 2 does seem 
to be approaching 0, as Dirichlet's Theorem predicts. The graph of f{x) with ^40 



IS 




Figure 4: Graph of f{x) and a sine series approximation of f{x). 
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2.2 Exercises in Fourier series using SAGE 

1. Let 

( 2, 3/2 < X < 3, 
f(x) = < 0, < X < 3/2, 
[ -1, -3<x < 0, 

period 6. 

• Find the Fourier series of f{x) and graph the function the FS 
converges to, — 3 < x < 3. 

• Write down the series in summation notation and compute the 
first 3 non-zero terms. 

solution: Of course the Fourier series of f{x) with period 2L is 



/ W ~ y + Z^K cos(-^) + bn sm(-^)J 

n=l 

where a„ and bn are 

1 TiTTX 

On = - y C0S(— ^)/(x) dx, 



1 ,nTcx 



^ , sin(——)f(x)dx. 
L J-L L ' 

Generally, the following functions compute the Fourier series coeffi- 
cients of /(x): 



def f ourier.series_cos_coef f (f cn,n,L) : # n>= 
lowerlimit = -L 
upperlimit = L 

an = maxima C ' tldef int C°/,s*cos C°/,s*n*x/°/,s) , x ,°/,s ,7,s) /L '"/, (fen, maxima (pi) ,L, lowerlimit , upperlimit) ) 
return str (an) .replace ("'/," , " " ) 

def f ourier_series_sin_coef f (f cn,n, L) : # n > 
lowerlimit = -L 
upperlimit = L 

bn = maxima C ' tldef int C°/,s*cos C°/,s*n*x/°/,s) , x ,-°/,s , °/,s) /L '/i (fen, maxima (pi) , L, lowerlimit , upper limit) ) 
return str (bn) .replace ("%" , " " ) 



However, Maxima's support for piecewise defined functions is rudimen- 
tary and the above functions will not give us what we want. So we 
compute them "by hand" but with some help from SAGE [S] (and 
Maxima, whichis included with SAGE): 
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## n > FS cosine coeff 

fen = 2; lowerlimit = 3/2; upperlimit =3; L = 3 

anl = maximaC 'tldef Int C7,s*cosC°/,s*n*x/7,s) ,x ,7,s ,°/,s) ' '/, (fen, maxima (pi) ,L, lowerlimit , upperlimit) ) ; anl 

# 6*sinC'/.pi*n)/(7,pi*n) - 6*sinC7opi*n/2)/C7.pi*n) 
fen = -1; lowerlimit = -3; upperlimit =0; L = 3 

an2 = maxima(*tldefintCXs*eos('/is*n*x/iis) ,x,%s ,'/,s) 'VtCf en,maximaCpi) ,L, lowerlimit, upperlimit)) ; an2 

# -3*sinC7.pi*n)/(7.pi*n) 
an = Canl+an2)/L 

# (3*sinC7.pi*n)/C7.pi*n) - 6*sinC7.pi*n/2)/Cy,pi*n))/3 



In other words, for n > 0, = 2 ^^^^^^^^^ 



## n - FS cosine coeff 

fen = 2; lowerlimit = 3/2; upperlimit =3; L = 3 

anl = maximaC'tldefintCXs*eosCHs*0*x/%s) ,x.,°/tB ,7,8) '7t(f en,maximaCpi) ,L, lowerlimit, upperlimit)) ; anl 

# 3 

fen = -1; lowerlimit = -3; upperlimit =0; L = 3 

an2 = maximaC'tldefintCXs*eos(Hs*0*x/%s) ,x.,°/tB ,7,8) '7,Cf en,maximaCpi) ,L, lowerlimit, upperlimit)) ; an2 

# -3 

an = Canl+an2)/L 

# 



In other words, ao = 0. 



## n > FS sine coeff 

fen = 2; lowerlimit = 3/2; upperlimit =3; L = 3 

bnl = maximaC 'tldef int CTis*sinC/is*n*x/5is) ,x,%s,y,s) 'y,Cf en,maximaCpi) ,L, lowerlimit, upperlimit)) ; bnl 

# 6*cos(7pi*n/2)/(7.pi*n) - 6*cosCXpi*n)/C7.pi*n) 
fen - -1; lowerlimit = -3; upperlimit =0; L = 3 

bn2 = maxima C ' tldef int C7,s*sinC7,s*n*x/7.s) , x ,7,s ,7,s) ' 7o( fen, maxima (pi) ,L, lowerlimit .upperlimit) ) ; bn2 

# 3/C7pi*n) - 3*cos(7,pi*n)/(7.pi*n) 
bn = (bnl+bn2)/L 

# ( - 9*cos(7.pi*n)/(7.pi*n) + 6*cos(7.pi*n/2)/(7.pi*n) + 3/(7.pi*n))/3 



In other words, bn = ^ +2cos(7rn/2) ^ 
The Fourier series is therefore 



N v~^r/r.sin(7rn/2) , .tittx. A — 3(— l)"" + 2cos(7rn/2), . .uttx 

- Ei(2^^)-«(^)+( )«-(^ 

n=l 



sage: Pi = RR(pi) 

sage: bn = lambda n: C - 9*cosCPi*n)/CPi*n) + 6*cosCPi*n/2)/CPi*n) + 3/CPi*n))/3 
sage: bn(l) ; bn(2) ; bn(3) 

1.2732395447351628 
-0 . 63661977236758127 
. 42441318157838753 



sage: an = lambda n: 2* (sin(Pi*n/2))/CPi*n) 
sage: anCl) ; anC2) ; anC3) 

0.63661977236768138 

. 00000000000000003898171832619376B 

-0.21220669078919379 
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Here are the first few numerical values of these coefficients: 



ai = ^ = 0.6366197723675813..., 

02 = 0, 

ag = -A = -0.2122065907891937..., 
hi = ^ = 1.273239544735162..., 

62 = V = -0.6366197723675812..., 

63 = 4 = 0.4244131815783875... . 



♦ 



Figure 5: Graph of Fourier series of f{x). 
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Let 

2, 3/2 < t < 3, 



0, 0<t< 3/2, 



• Find the Cosine series of f{x) (period 6) and graph the function 
the CS converges to, — 3 < x < 3. 

• Write down the series in summation notation and compute the 
first 3 non-zero terms. 

solution: Of course the Cosine series of f{x) with period 2L is 



fix) ~ Y + 2^anCos(-^), 

where a„ is 



2 

n=l 



J cos(^^)/(x) dx. 



A simple Python program: 



def cosine_series_coeff Cfcn,ii,L) : # ii>= 
lowerlimit = 

upperlimit = L 

an = maxima C ' 2*tldef int C°/oS*cos(lis*n*x/%s) ,x,iis ,y.s)/L *X (fen, maxima Cpi) , L, lowerlimit , upperlimit ) ) 
return strCan) .replace ("%" , " ") 



It was noted above that this program will not help us, for the type of 
function wc arc dealing with here. So, we have SAGE do the compu- 
tations "piece-by- piece" : 

## n > 

fen = 2; lowerlimit = 3/2; upperlimit =3; L = 3 

anl = maximaC 'tldef int C*;(s*cos(*;(s*n*x/%s) ,^,7,5 ,'/ts') 'VtCf cn,maximaCpi) ,L, lowerlimit , upperlimit) ) ; anl 

# 6*sinC%pi*n)/(7.pi*n) - 6*sinC%pi*n/2) /C%pi*n) 
an = C2/L)*anl 

an 

# 2*(6*sin(7.pi*n)/("/,pi*n) - 6*sin(7.pi*n/2) / C7.pi*n) )/3 



In other words, a„ = 4£Ei![!iM Xo find the 0-th coefficient, use the 

' "■ Trn ' 

commands 



## n - 

anl = maximaC 'tldef int Cyis*cos(yis*0*x/7.s) ,x,7s,7,s) ' 7tCf cn,maximaCpi) ,L, lowerlimit , upperlimit) ) ; anl 

# 3 

an = (2/L)*anl 

an 

# 2 

# aO - 2 



28 



or the command 



sage 
sage 
sage 
sage 



fl = lambda x:2 
f2 = lambda x:0 



f = PiecewiseC[[(0,3/2),f2],[(3/2,3),fl]]) 
f . cosine_series_coef f icient (0,3) 



In other words, ao — 2. 



The Cosine series is therefore 



oo 



sin(7rn/2) 




n=l 



sage: an = lambda n:4*sin(Pi*n/2)/(Pi *n) 
sage: an(l) ; an(2) ; an(3) 

1.2732395447351628 

. 000000000000000077963436650387510 

-0.42441318157838759 

So, ai = ^ = 1.273239544735162..., 02 = 0, 03 = -f . 
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Figure 6: Graph of Cosine series of /(x). 



Let 



fix) 



2, 3/2<t<3, 
0, < t < 3/2, 



Find the Sine series of f{x) (period 6) and graph the function the 
FS converges to, — 6 < x < 6. 

Write down the series in summation notation and compute the 
first 3 non-zero terms. 



We leave this one as an exercise for the reader! 
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2.3 Application to the heat equation 

The heat equation with zero ends boundary conditions models the tempera- 
ture of an (insulated) wire of length L: 

j^d-^u(x,t) du{x,t) 

u{0,t) ^u{L,t) = 0. 

Here u{x,t) denotes the temperature at a point x on the wire at time t. 
The initial temperature f{x) is specified by the equation 

u{x,0) = f{x). 

Method: 

• Find the sine series of f{x): 

oo 

f{x)r.J2bn{f)sini^), 

n=l 

• The solution is 



n=l 

Example 14 Let 
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f(^, _ / -1' < X < 3/2, 
' \ 2, 3/2 <x< 3. 

Then L = 3 and 

2 /"^ 

K{f) = I /(x) sin(n7rx/3) (ix. 
-3 Jo 

Thus 



/(x) ~ 6i(/)sin(x7r/3) + 62(/)sin(2x7r/3) + ... = - sin(7rx/3)) - - sin(27rx/3) + ... . 

vr vr 

The function /(x), and some of the partial sums of its sine series, looks like 
Figure 




Figure 7: /(x) and two sine series approximations, 6*10, 5*30. 
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As you can see, taking more and more terms gives functions which better and 
better approximate f{x). 



sage: R = PolynomialRingCQQ, "x") ; x = R.genO 

sage: fl = -x~0; ±2 = 2*x~0 

sage: f = Piecewise ([ [(-3,-3/2) ,-f 2] , [(-3/2,0) ,-fl] , [(0 ,3/2) ,fl] , [(3/2,3) ,f 2] ] ) 

sage: fslO = f . f ourier_series_partial_sum(10 ,3) 

sage: FSIO = lambda t :RR(sage_eval (f slO .replaceC'x" ,str (t) ) ) ) 

sage: PfslO = plot (FSIO, -3, 3) 

sage: Pf = f.plotO 

sage: show(Pf +Pf slO) 

sage: fs30 = f .f ourier_series_partial_sum(30,3) 

sage: FS30 = lambda t : RR(sage_eval (fs30 .replaceC'x" ,str(t) )) ) 

sage: Pfs30 = plot(FS30,-3,3) 

sage: showCPf +Pf slO+Pf s30) 

The solution to the heat equation, therefore, is 

oo 

u{x,t) = J2 KU) sin(^) eM-kC^?t)- 



The heat equation with insulated ends boundary conditions models the 
temperature of an (insulated) wire of length L: 

j^ d'^u(x,t) du(x,t) 

u.,{0,t) =u,,{L,t) = 0. 

Here Ux{x, t) denotes the partial derivative of the temperature at a point x on 
the wire at time t. The initial temperature f{x) is specified by the equation 
u{x,0)^f{x). 
Method: 

• Find the cosine series of f{x): 

oo 

fix) ^ Y + X] ^^^i-f) 

n=l 

• The solution is 

oo 

^) y + 2^ «n(/) cos(-^)) exp{-k{—) t). 

n=l 
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Example 15 Let 



Then L = tt and 



for n > and oq = 1. 
Thus 



fix) 



-1, 0<x<3/2, 
2, 7t/2<x < 3. 



2 

an(/) = ^ / f{x)cos{mTx/3)dx, 

-J JO 



/(x) cos(x/3)+a2(/) cos(2z/3)+...l/2 cos(7rx/3)H — cos(37rx/3)+.. 

2 TT TT 

The piecewise constant function f{x), and some of the partial sums of its cosine 
series, looks like Figure |HI 




Figure 8: f{x) and two cosine series approximations. 
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sage: R = PolynomialRingCQQ, "x") ; x = R.genO 

sage: fl = -x~0; ±2 = 2*x~0 

sage: f = PiecewiseC [[(-3,-3/2) ,f2] , [(-3/2,3/2) ,fl] , [(3/2,3) ,f2]] ) 

sage: fslO = f . f ourier_series_partial_sum(10 ,3) 

sage: FSIO = lambda t :RR(sage_eval (f slO .replaceC'x" , str (t) ) ) ) 

sage: fs30 = f .f ourier_series_partial_sum(30,3) 

sage: FS30 = lambda t :RR(sage_eval (fs30 .replaceC'x" , str (t) )) ) 

sage: Pf = f.plotO 

sage: Pfs30 = plot(FS30,-3,3) 

sage: PfslO = plot(FS10,-3,3) 

sage: show(Pf +Pf slO+Pf s30) 

As you can see, taking more and more terms gives functions which better and 
better approximate f{x). 

The solution to the heat equation, therefore, is 

oo 

■"(^^ *) = y + XI «"(/) cos(^^) exp(-fe(^)2t). 

n=l 



Explanation: 

Where docs this solution come from? It comes from the method of sepa- 
ration of variables and the superposution principle. Here is a short explana- 
tion. We shall only discuss the "zero ends" case (the "insulated ends" case 
is similar). 

First, assume the solution to the PDE k^^^^ = has the "factored" 

form 



u{x,t) = X{x)T{t), 

for some (unknown) functions X, T. If this function solves the PDE then it 
must satisfy kX"{x)T{t) = X{x)T'{t), or 

X"{x) _ 1 T'{t) 
X{x) ~ k T{t) ' 

Since x, t are independent variables, these quotients must be constant. In 
other words, there must be a constant C such that 

T'(f) 

J^ = kC, X"{x)-CX{x) = 0. 
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Now we have reduced the problem of solving the one PDE to two ODEs 
(which is good), but with the price that we have introduced a constant which 
we don't know, namely C (which maybe isn't so good). The first ODE is 
easy to solve: 

T(t) = Aie''^\ 

for some constant Ai. To obtain physically meaningful solutions, we do not 
want the temperature of the wire to become unbounded as time increased 
(otherwise, the wire would simply melt eventually). Therefore, we may as- 
sume here that C < 0. It is best to analyse two cases now: 

Case C = 0: This implies X{x) = A2 + A^x, for some constants A2, A3. 
Therefore 

u{x, t) = Ai{A2 + A3X) = y + box, 

where (for reasons explained later) A1A2 has been renamed ^ ^md A1A3 has 
been renamed 60 • 

Case C < 0: Write (for convenience) C = — r^, for some r > 0. The ODE 
for X implies X{x) = ^2 cos(ra;) + y43sin(rx), for some constants ^2, A3. 
Therefore 



u{x,t) = Aic ^"^ *(A2Cos(rx) +y43sin(rx)) = (acos(rx) + 6sin(rx))e ^"^ *, 

where A1A2 has been renamed a and AiA^ has been renamed h. 

These are the solutions of the heat equation which can be written in 
factored form. By superposition, "the general solution" is a sum of these: 

u{x,t) = |i + 60a; + X]r=i(«nCos(r„x) + 6„sin(r„x))e"'"'"* 

= Y + + (oi cos(rix) + 61 sin(rix))e~'^''i* (7) 
+ (02 cos(r2x) + 62 sin(r2x))e~^^2* + 

for some Cj, bi, rj. We may order the r^'s to be strictly increasing if we like. 

We have not yet used the IC u{x, 0) = f{x) or the BCs m(0, t) = u{L, t) = 
0. We do that next. 

What do the BCs tell us? Plugging in x = into (jT)) gives 
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= u{0, t) = ^ + J2 «ne-'^'^"* = ^ + a,e-''"' + a^e-'"-"' + ... . 

n=l 

These exponential functions are linearly independent, so ao = 0, ai = 0, 
a2 = 0, ... . This implies 

u{x,t) = 6oa;+^6nSin(r„x)e"^''"* = box+bi sin(rix)e"'"'i*+62 sin(r2x)e"'"'2*+... . 

n=l 

Plugging in X = L into this gives 

= u{L, t) = boL + '^ bn sin(r„L)e"^''"*. 

n=l 

Again, exponential functions are linearly independent, so bo = 0, 6„ sin(r„L) 
for n = 1, 2, .... In other to get a non-trivial solution to the PDE, we don't 
want bn = 0, so sin(r„L) = 0. This forces r„L to be a multiple of vr, say 
r„ = rnr/L. This gives 



= V'6„sin(^x)e ''^t? = 6i sin(Yx))e ^^^^'*+62 sin(^a;))e ''(^)'*+. 

1j 1j 

n=l 

(8) 

for some 6j's. This was discovered by Fourier. 

There is one remaining condition which our solution u{x,t) must satisfy. 
What does the IC tell us? Plugging t = into (jHI) gives 



J[x) = u[x,0) = 2_^bnSm{—x) = Oism(— x)) +02sm(— xj) + ... . 



n=l 



In other words, if f{x) is given as a sum of these sine functions, or if we can 
somehow express f{x) as a sum of sine functions, then we can solve the heat 
equation. In fact there is a formula for these coefficients bn (which Fourier 
did not know at the time): 



bn = Y [ f{x) cos(^x)dx. 
Jo ^ 



It is this formula which is used in the solutions above. 
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Example 16 Let = 1, let 

m = 

and let g{x) = 0. Then L = 3 and 

Mf) =i/o/(x)sin(^)d: 



-1, 0<a;<3/2, 
2, 3/2<x<3. 



-2 

nn 



2 cos(n7r)— 3 cos(l/2 n7r)+l 

Thus 



fix) ~ 6i(/) sin(x7r/3) + hif) sin(2x7r/3) + ... 
= f sin(^x/3)) - f sin(27rx/3) + ... . 

The function f{x), and some of the partial sums of its sine series, looks like Figure 
[7| The solution is 



= I sin(7ra;/3))e-(^/3)'* - f sin(27ra;/3)e-(2V3)2t + 
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Figure 9: Plot of /(x) and u{x, t), for t = 0.0, 0.3, 0.6, 0.9. The first plot uses 
the partial sums 6*50 of the FS for /(x); the second plot uses the Cesaro- 
filtered partial sums 5*^ of the FS for f{x). 



The following SAGE code for the plot above is very time-consuming: 



sage: fl = lambda x:-2 

sage: f2 = lambda x:l 

sage: f3 = lambda x:-l 

sage: f4 = lambda x:2 
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sage: f = Piecewise ([ [(-3,-3/2) ,fl] , [(-3/2,0) ,f 2] , [(0,3/2) ,f 3] , [(3/2,3) ,f 4] ] ) 

sage: N = 50 

sage: t = 0.0; F = [RR(exp(-(j*pi/3) "2*t) ) for j in range(N)] 

sage: PC = f .plot_f ourier_series_partial_suin_f iltered(N,pi,F,-4,4) 

sage: t = 0.3; F = [RR(exp(-(j*pi/3)~2*t)) for j in range(N)] 

sage: PI = f .plot_fourier_series_partial_sum_filtered(N,pi,F,-4,4) 

sage: t = 0.6; F = [RR(exp(-(j*pi/3)~2*t)) for j in range(N)] 

sage: P2 = f .plot_f ourier_series_partial_suin_f iltered(N,pi,F,-4,4) 

sage: t = 0.9; F = [RR(exp(-(j*pi/3) "2*t) ) for j in range(N)] 

sage: P3 = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: P = f .plot (rgbcolor=(0 . 8,0 . 1 ,0 .3) , plot_points=40) 

sage: show(P+P0+Pl+P2) 

sage : N = 50 

sage: t = 0.0; F = [RR((l-j/M)*exp(-(j*pi/3)~2*t)) for j in range(N)] 

sage: PcO = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: t = 0.3; F = [RR( (1-j/N) *exp(-( j *pi/3) '2*t) ) for j in range(N)] 

sage: Pel = f .plot_f ourier_series_partial_suin_f iltered(N,pi,F,-4,4) 

sage: t = 0.6; F = [RR( (1-j/N) *exp(-( j *pi/3) '2*t) ) for j in range(N)] 

sage: Pc2 = f .plot_f ourier_series_partial_suin_f iltered(N,pi,F,-4,4) 

sage: t = 0.9; F = [RR((l-j/N)*exp(-(j*pi/3)"2*t)) for j in raiige(N)] 

sage: Pc3 = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: show(P+Pc0+Pcl+Pc2) 



2.4 Application to Schodinger's equation 

The one-dimensional Schrodinger equation for a free particle is 

.j^d'^ip{x,t) dilj{x,t) 
dx^ dt 

where /c > is a constant (involving Planck's constant and the mass of 
the particle) and i — y/—! as usual. The solution ip is called the wave 
function describing instantaneous "state" of the particle. For the analog 
in 3 dimensions (which is the one actually used by physicists - the one- 
dimensional version we arc dealing with is a simplified mathematical model), 
one can interpret the square of the absolute value of the wave function as the 
probability density function for the particle to be found at a point in space. 
In other words, {x,t)f dx is the probability of finding the particle in the 
"volume dx" surrounding the position x, at time t. 

If we restrict the particle to a "box" then (for our simplied one-dimensional 
quantum-mechanical model) we can impose a boundary condition of the form 
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ij{0,t)^iP{L,t) = 0, 
and an initial condition of the form 



il^lx, 0) = f{x), < X < L. 

Here / is a function (sometimes simply denoted iIj{x)) which is normalized 
so that 

r\f{x)\'dx^i. 

Jo 

If \ilj{x,t)f represents a pdf of finding a particle "at x" at time t then 
/o \f{^)\'^dx represents the probability of finding the particle somewhere 
in the "box" initially, which is of course 1. 

Method: 

• Find the sine series of f{x): 

oo 

/(x)~5]M/)sm(^), 

n=l 

• The solution is 



n=l 

Each of the terms 



ipn\x,t) = f)„sm(-^)exp(-zA;(— ) t). 

is called a standing wave (though in this case sometimes 6„ is chosen so that 
j^\^^{x,t)\''dx^l). 



Example 17 Let 



/(X) = { - 



1, 0<x<l/2, 
1/2 < X < 1. 
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Then L = 1 and 



, / 2 /"-^ „. , . .mix. , 1 / -, „ t'^T^^ , XN 

= T / /WSin(— — jdx = — (-1 + 2cos(— ) - cos(n7r)). 
1 7o 1 nTT 2 



Thus 



fix) ~ sin(7rx) + sin(27rx) + ... 

= En + 2cos(^) - cos(n7r)) • sin(n7rx) . 

Taking more and more terms gives functions which better and better approxi- 
mate j{x). The solution to Schrodinger's equation, therefore, is 

OO J 

ip{x,t) = ^ — (— 1 + 2cos(— ) — cos(n7r)) ■ sin(n7ra;) • ex.p{—ik{mT)^t). 

n=l 



Explanation: 

Where does this solution come from? It comes from the method of separa- 
tion of variables and the superposution principle. Here is a short explanation. 

First, assume the solution to the PDE ik^^^ = has the "fac- 

tored" form 

ij{x,t)^X{x)T{t), 

for some (unknown) functions X, T. If this function solves the PDE then it 
must satisfy kX"{x)T{t) = X{x)T'{t), or 

X"{x) _ 1 T'{t) 
X(x) ~ ik T(t) ' 

Since x, t are independent variables, these quotients must be constant. In 
other words, there must be a constant C such that 

T'(f) 

-Y^ = ^kC, X"{x)-CX{x)=0. 

Now we have reduced the problem of solving the one PDE to two ODEs 

(which is good), but with the price that we have introduced a constant which 
wc don't know, namely C (which maybe isn't so good). The first ODE is 
easy to solve: 
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T(t) = Aie'^^\ 

for some constant Ai. It remains to "determine" C. 

Case C > 0: Write (for convenience) C = r^, for some r > 0. The ODE 
for X imphes X{x) = A2 exp(ra;) + A3 exp(— rx), for some constants A2, A^. 
Therefore 



ipix^t) = Axe *(y42exp(ra;)+y43exp(-ra;)) = (aexp(rx)+6exp(-ra;))e 

where ^41^42 has been renamed a and Ai^a has been renamed h. This will 
not match the boundary conditions unless a and h are both 0. 

Case C = 0: This implies X[x) = A^^ A-^x^ for some constants A2^ A3. 
Therefore 

t) = A\{A2 + A3X) = a + 5x, 

where AxA^ has been renamed a and ^1^3 has been renamed h. This will 
not match the boundary conditions unless a and h are both 0. 

Case C < 0: Write (for convenience) C = — r^, for some r > 0. The ODE 
for X implies X{x) = ^2 cos(rx) + A3sin(rx), for some constants A2, A3. 
Therefore 



ip{x,t) = Aie *(A2Cos(rx) + A3sin(rx)) = (acos(rx) + 6sin(rx))e *, 

where A1A2 has been renamed a and A1A3 has been renamed b. This will 
not match the boundary conditions unless a = and r = ^ 

These are the solutions of the heat equation which can be written in 
factored form. By superposition, "the general solution" is a sum of these: 

^(2^5 1) = Z]r=i(an cos(r„a;) + 6„ sin(r„x))e"*''''"* 
= 61 sin(ria;)e~*'^'''i* + 62 sin(r2a;)e~*'^^'2* + 

for some 6„, where = Note the similarity with Fourier's solution to 
the heat equation. 

There is one remaining condition which our solution il){x,t) must satisfy. 
We have not yet used the IC tfj^x, 0) = f{x). We do that next. 
Plugging t = into gives 
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oo 

, . ,n7r ^ , . . TT , , , . ,271 



f{x) = i>{x, 0) = 2_^bn sm{—x) = bi sin(— a;)) + 62 sm{—x)) + ... . 

In other words, if f{x) is given as a sum of these sine functions, or if we 
can somehow express f{x) as a sum of sine functions, then we can solve 
Schrodinger's equation. In fact there is a formula for these coefficients 

2 '•^ 



/ f{x) cosC^x)dx. 
Jo L 



" L 

It is this formula which is used in the solutions above. 



2.5 Application to the wave equation 

The wave equation with zero ends boundary conditions models the motion 
of a (perfectly elastic) guitar string of length L: 



2 d w{x, 

w{^,t) ^ w{L,t) =0. 



Here w(a;, t) denotes the displacement from rest of a point x on the string at 
time t. The initial displacement f{x) and initial velocity g{x) are specified 
by the equations 

w(a;,0) = /(x), Wt{x,Gi) = g{x). 

Method: 

• Find the sine series of f{x) and g{x): 

00 mix °° TT 

fi.x) ~ Kif) sin(-^), g{x) ~ sm(-^). 

n=l n=l 

• The solution is 

w{x,t) = V 6n / COS — — + ^i^sm — — sm -- . 
^-^ L nna L L 

n=l 
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A special case: When there is no initial velocity then g = and the 
solution to the wave equation, therefore, is 



n=l 



-1, 0<x<7r/2, 
2, 7r/2 < X < vr. 



Example 18 Let a = 1, let 

fix)-- 

and let g(x) = 0. Then L = vr, bn{g) = 0, and 

bnif) = f /J" /(x) sm{nx)dx 

2 2 cos(n7r)— 3 cos{l/2 n7r)+l 



nn 



Thus 



/(x) - sm(x) + 62(/) sm(2x) + ... 



= I sm(a;) - | sin(2a;) + J;: sin(3a;) + .... 

The function f{x), and some of the partial sums of its sine series, looks like Figure 
[7| The solution is 



w{x, t) = Yl'^^i hn{f) cos(nt) sin(nx) 

= ;| cos(t) sin(x) — ^ cos(2t) sin(22;) + ^ cos(3t) sin(3x) + ... . 
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Figure 10: Plot of /(x) and w{x,t), for t = 0.0,0.3,0.6,0.9. The first plot 
uses the partial sums S25 of the FS for /(x); the second plot uses the Cesaro- 
filtered partial sums 5*^ of the FS for f{x). 

The following SAGE code for the plot above is very time-consuming: 



sage: fl = lambda x:-2 

sage: f2 = lambda x:l 

sage: f3 = lambda x:-l 

sage: f4 = lambda x:2 
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sage: f = Piecewise ( [ [(-pi ,-pi/2) ,f 1] , [(-pi/2,0) ,f 2] , [(0,pi/2) ,f 3] , [(pi/2,pi) ,f 4] ] ) 

sage : N = 25 

sage: t = 0.0; F = [RR(cos((j+l)*t)) for j in range(N)] 

sage: PC = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: t = 0.3; F = [RR(cos((j+l)*t)) for j in range(N)] 

sage: PI = f .plot_fourier_series_partial_sum_filtered(N,pi,F,-4,4) 

sage: t = 0.6; F = [RR(cos ( ( j+1) *t) ) for j in range(N)] 

sage: P2 = f .plot_f ourier_series_partial_suin_f iltered(N,pi,F,-4,4) 

sage: t = 0.9; F = [RR(cos((j+l)*t)) for j in range(N)] 

sage: P3 = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: P = f .plot (rgbcolor=(0 . 8,0 . 1 ,0 .3) , plot_points=40) 

sage: show(P+P0+Pl+P2) 

sage: fl = lambda x:-2 

sage: f2 = lambda x:l 

sage: f3 = lambda x:-l 

sage: f4 = lambda x:2 

sage: f = Piecewise ([ [(-pi , -pi/2) ,fl] , [(-pi/2,0) ,f 2] , [(0, pi/2) ,f 3] , [(pi/2, pi) ,f 4] ] ) 

sage : N = 25 

sage: t = 0.0; F = [RR((l-j/N)*cos((j+l)*t)) for j in range(N)] 

sage: PO = f .plot_fourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: t = 0.3; F = [RR((l-j/N)*cos((j+l)*t)) for j in range(N)] 

sage: PI = f .plot_f ourier_series_partial_sum_f iltered(N,pi,F,-4,4) 

sage: t = 0.6; F = [RR((l-j/N)*cos((j+l)*t)) for j in range(N)] 

sage: P2 = f .plot_f ourier_series_partial_sum_f iltered(N,pi ,F, -4,4) 

sage: t = 0.9; F = [RR((l-j/N)*cos((j+l)*t)) for j in range(N)] 

sage: P3 = f .plot_fourier_series_partial_sum_filtered(N,pi,F,-4,4) 

sage: P = f .plot (rgbcolor=(0 .8,0 . 1 ,0 .3) , plot_points=40) 

sage: show(P+P0+Pl+P2) 



3 The Discrete Fourier transform 

Let us first "discretize" the integral for the k-th coefficient of the Fourier 
series and use that as a basis for defining the DFT. Using the "left-hand 
Riemann sum" approximation for the integral using N subdivisions, we have 

^ ^ f{Pj/N)e-^^'''^'^m (10) 
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This motivates the following definition. 



Definition 19 The N -point discrete Fourier transform (or DFT) of the vector 
/ = (/o,...,/iv-i)eC^ is 



7V-1 N-1 

DFT{f)k = fk=Yl /,e-2-*^^/^ = ^ fjW''', 

j=0 j=0 

where W = e^^^l^ . 

The normalized N -point discrete Fourier transform (or NDFT) of the vector 
/= (/o,-,/iv-i)eC^ is 



NDFT{f)k 



1 



N-l 

Hi 

j=0 



N-l 



-2irikj/N 



-1= y f^w^r 



Note that the powers of W are equally distributed points on the unit 
circle. 

Both the DFT and NDFT define linear transformations and 

— * 

therefore can be described by matrices. If we regard the vector / as a column 
vector then the matrix for the DFT is: 





( 1 


1 




1 


W 




1 






I 1 





■2(Ar-i) 



W 



_(jV-l)(iV-l) 



Note that this is a symmetric matrix. Similarly, the matrix for the NDFT is: 



/ 1 


1 


1 


w 


1 




\1 





1 \ 

■iV-1 



W 



w 



iN-l){N-l) 



J 



Example 20 Let N ^10. The DFT of 

f = (1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10) e c 



10 
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Fio/ = (1,0,0,0,0,0,0,0,0,0). The DFT of 



f= (1/10, 0, 0, 0, 0, 0, 0, 0, 0, 0) e 0^ 

IS Fwf = (1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10). 



This was computed using the SAGE commands 

sage: J = range (10) 
sage: A = [1/10 for j in J] 
sage: s = IndexedSequence(A, J) 
sage: s.dftO 

Indexed sequence: [1, 0, 0, 0, 0, 0, 0, 0, 0, 0] 
indexed by [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] 

There is an analog of the orthogonality used in Theorem 0] above. 
Lemma 21 We have 



Before proving this algebraically, I claim that this is "geometrically ob- 
vious." To see this, recall that the average of any points in the plane - 
whether written as vectors or as complex numbers - is simply the "center 
of gravity" of points, regarded as equally weighted point masses. The sum 
above is (if does not divide k) the "center of gravity" of a collection of 
point masses which are equi-distributed about the unit circle. This center of 
gravity must be the origin. On the other hand, if N\k then all the points are 
concentrated at 1, so the total mass is in that case. 

proof: If W 7^ 1 then we have 



IfW^ = 1 then we "£^=0 = Note w'' = I if and only if N\k. □ 

As a corollary of this lemma, we see that the complex matrix F/v is 
"orthogonal" (this is a technical term we need to define). A real square 




j=0 



otherwise. 
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matrix is called orthogonal if row i is orthogonal to column j for all i ^ j. 
Here when we say two real vectors are orthogonal of course we mean that 
they are non-zero vectors and that their dot product is 0. The definition 
for complex matrices is a bit different. We first define on the Hermitian 
inner product (or just inner product, for short): 

N-l 
j=0 

We say two complex vectors are orthogonal if they are non-zero and their 
inner product is zero. A complex square matrix is called orthogonal if row i 
is orthogonal to column j for all i ^ j- 

Lemma 22 Fjv is orthogonal. 

proof: The k-t\i row of Fn is the vector (W )j=o,...,N-i-, and the com- 
plex conjugate of this vector is the vector iyV''^~^^-')j=Q^...,N-i = {W ^^'')j=o,...,N-i, 
so 

iV-l 

((row k of Fn), (row i of Fn)) = ^ g, 

j=0 

provided W 7^ 1, which is true if and only if does not divide k — i. d 
Note that this matrix F/v is not "real orthogonal" : 

N-l 

(row k of Fn) ■ (row i of Fn) = ^ = q, 

j=0 

if and only if does not divide k + i — 2, by Lemma 1^ 

Here's another matrix calculation based on Lemma El If Cfc denotes 
the standard basis vector of whose fc-th coordinate is 1 and all other 
coordinates are 0, then 

DFT{ek) = (Fw)fe = k*'' column of F^, 

so 
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DFT{DFT{ek)) = DFT{k^'^ column of F, 



N 



( 



(1** row of Fiv 
(2"°' row of 



AT 



■ (A;*^ column of Fat) \ 

■ {k^^ column of F/v) 



y {N^^ row of Fat) • (A;*^ column of Fat) y 



/ (1^* row of F^) 
(2"'^ row of Fat) 



(fc*'^ row of Fiv 
(A;*'^ row of F, 



\^ (A^*'^ row of F^) ■ (A;*'' row of F^) j 

since the matrix Fjy is symmetric. The "almost orthogonality of the rows of 
Fat" discussed above imphes that this last vector of dot products is = Nc-ki 
where by —k in the subscript of we mean the representative of the residue 
class of —k (mod A^) in the interval < —k < N — 1. 

The motivation behind the definition of the normalized DFT is the fol- 
lowing computation: 



NDFT{NDFT{f))u = ^ Ef=o' NDFTif),W 



1_ 

N 



f- 



EN-l f 
£=0 nl^ 



N-l 
j=0 



W 



{k+e)j 



'111 



where of course by —k in the subscript of f^k we mean the representative 
of the residue class of —A; (mod A^) in the interval < —k < N — 1. 
To be precise, if neg : denotes the negation operator, send- 

ing (/o,/i,...,/iv-i) to (/_o,/-i,...,/i-iv) then NDFT^ = neg. Note that 
neg flips the last — 1 coordinates of / about their midpoint. For ex- 
ample, if A^ = 5 then neg(l, 2, 3, 4, 5) = (1,5,4,3,2) and if A^ = 6 then 
neg(2, 1,3, -1,4, 7) = (2, 7, 4, -1,3,1). 

Because of the computation (jUj), it follows that 



NDFT~\f)k = NDFT{f). 



~kj 



j=0 



N-l 



(12) 



3=0 
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or NDFT~^ = neg o NDFT = NDFT o neg. Likewise, 



^ N-l N-1 
j=0 j=0 

or DFT-^ = N-^ ■ neg o DFT = N'^ ■ DFT o neg. 
Example 23 Let N = 4, so 

/II 1 1 \ 




\ 1 i -I -i J 
and let f = (1, 0, 0, 1). We compute 



NDFTif) = ±Fj= i 



/II 1 

1 -i -1 

1 -1 1 

\1 i -1 



1 \ 



/ 1 \ 





VI/ 



/ 1 \ 

(1 + 0/2 


V (1-0/2 / 



Call this latter vector g. We compute 



NDFT o neg( 



/II 1 

1 -I -1 

1 -1 1 

\l i -I 



1 \ 



/ 



(1 + 0/2 


V (1-0/2 / 



/2\ 




V2/ 



as desired. 



3.1 Eigenvalues and eigenvectors of the DFT 

We shall take an aside to try to address the problem of computing eigenvalues 
and eigenvectors of DFT = F^, at least in a special case. This shall not be 
used in other parts of the course - think of this as "for your cultural benefit" . 

It also follows from this computation (fTT|) that NDFT^ acts as the iden- 
tity. If A is any square matrix satisfying A"^ = I then any eigenvalue A 
of A must be an m-th root of unity (indeed, Av = Xv, for some non-zero 
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eigenvector v, so v = A"^v = A"* ^Av = = \A"^ = ... = X"^v, 

so A"^ = 1). Due to this fact, it follows that the only possible eigenvalues of 
NDFT are 1,-1, 

It is not hard to describe the eigenspaces intrinsically. Let 

V = {f : Z/NZ C}, 
so y ^ via / ^ (/(O), /(I), fiN - 1)). Let 

Veven = {f & V \ f{-k) = f (k) , G Z/NZ} 

denote the subspace of even functions and 

Vodd = {fev\ f{-k) = -f{k), \Jk e Z/NZ} 

the subspace of odd functions. Note that the restriction of NDFT to Veven 
is order 2: for all / G Vgyen, we have NDFT^{f) = f. Likewise, for all 
/ G Vodd, we have NDFT\f) = -/. Let 

El = {NDFTif) + / I / G Veven}, 
E_i = {NDFTif) -fife I4.en}, 

E., = {iNDFTU) + f\fe Vodd}, 

E, = {iNDFTU) -fife Kdd}. 

In this notation, for each A G {±l,±i}, E\is the eigenspace of Gat having 
eigenvalue A. Moreover, according to Gooqj the columns of the matrix 

Mx = I + XGn + X'^Gl + A^G^ (14) 

form a basis for Ex. 

Note that the eigenvalues of the DFT = F/v are not the same as those of 
NDFT = Gn since G^ — "^Fn- The eigenvalues of the DFT must belong 

to y/N, -y/N, iy/N, -iy/N. 

Recall that for the Fourier transform on M the number A = 1 was an 
eigenvalue. The calculation there cannot be "discretized" since we had to 

^Good's definition of the NDFT is slightly different than ours. Essentially, where we 
have a weighted sum over powers of W, he has a weighted sum over powers of W. That 
changes the matrix AI\ slightly, so the one above is correct for us I think. 
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use the Residue Theorem from complex analysis. We have to take another 
approach. 

Let Tj/NX = {0, 1, 2, A^— 1}. This is a set but if we perform arithmetic 
(such as addition and multiplication) mod N, then this can regarded as an 
abelian group. 

First, assume that there is a function 

i : Z/NZ Z/NZ (15) 

with the property that i{jk) = i{j)i{k), for all j ^ and k ^ 0, and i{0) = 0. 
Let us also assume that the set 

kZ/NZ = {jk (mod N) \j^ 0, 1, - 1} 
which is a subset of Z/NZ, is the same as the set Z/NZ: 

kZ/NZ = Z/NZ, 0<k<N. (16) 

We will worry about whether this function exists or not and whether this 
set-theoretic property of true or not later. For now, let's compute the first 
component of it's DFT: 

N-l 

FDT(/)i = ^^(j)e-2"^/^ 

j=0 

where /= {i{0),i{l), ...,i{N - 1)) = (0, £(1), £(iV - 1)). Now make the 
substitution j = j'k, for some k ^ 0. 
We have 

FDT{e)^ = EJTo' ^OOe"'"^'/"^ 
= i{k)DFT{i)k. 

Putting these together, we get 

DFT{e) = DFT{i)i-l 

— * — * 

In other words, i. is an eigenvector with eigenvalue DFT{£)i. 
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Example 24 In general, let Cn denote a primitive N-th root of unity, for 
example Cn = e^^'^*/^ = W. The smallest field containing the rationals, Q, 
and this N-th root of unity is called the cyclotomic field, and is commonly 
denoted by Q,{Cn)- a set, Q(CAr) is simply the set of "polynomials in (j^ 
with coefficients in Q of degree < N — 1". 

Let N = 5 and i = (0, 1, —1, —1, 1). It can be checked that this defines a 
function as in 1115]} . 

We first check explicitly that assumption / Ii6|) is true: 

Z/5Z = [0,1,2,3,4] 
2Z/5Z = [0,2,4,1,3] 
3Z/5Z = [0,3,1,4,2] 
4Z/5Z = [0,4,3,2,1] 

Let 



where C5 



-CI 



/ 1 1 



1 \ 



1 Cs Cl Cl Cl 

1 cl cl Cs cl 

1 C5 Cs C5 C5 

*3 a2 



V 1 Cl Cl Cl Cs / 

Cs — 1 . The characteristic polynomial of this matrix is 



P5{x) = {x- V5){x + V5)^(x^ + 5). 

The roots of this polynomial are of course the eigenvalues: ±-\/5, ±i-\/5- Of 
course, the matrix F5 has eigenvectors with vector components in C. We 
shall try to express their components, if we can, algebraically in terms of 
the powers of the Cs- This is not a matter of necessity for us, but it can 
be convenient for doing certain calculations. For example, you don't have to 
worry about round- off errors messing up a calculation if you have an algebraic 
expression. 

It turns out that a/5 = -2C| - 2C| - 1, and -^5 = 2C| + 2C| + 1, can 
both be written as a linear combination of powers of (5, so it may come as 
no surprise that the components of their eigenvectors also can be written as 
a linear combination of powers of (5. In other words, if the eigenvalue is in 
QlCs) then one might expect to find eigenvector with components in Q(C5)- 

On the other hand, ±i\/5 do not belong to Q(C5)- So, we should expect 
that the eigenvectors in this case have components lying in some extension of 
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in Q(C5)- It turns out, all the eigenvectors have components in Q(C2o); the 
field extension of Q generated by the 20-th roots of unity. 
The eigenspace of Xi = is 2- dimensional, with basis 

(1, 0, -cl - C5 - 1, -Cl - C5 - 1, 0), (0, 1, -1, -1, 1), 

where —Q — — 1 = 0.6180... . The eigenspace of X2 = —V5 is 1- 
dimensional, with basis 



1^3 , K2 Ks 



1.2 1.3 



1 

2' 



where ICf + ld = —0.8090... . The eigenspace 0/A3 = i^/E is 1- dimensional, 
with basis 

(0, 1, C20 + C20~C20"~C20 + C20~2C20-1, "Cjo ~ C20 + C20 + C20 ~" Clo + ^(20 + 1, "1), 

where (lo + Clo " Clo - C20 + Clo " 2(20 - 1 = -3.520... . The eigenspace of 
A4 = i\/5 is 1- dimensional, with basis 

(0, 1, -C20 + C20 + C20~C20~C20 + 2C20-1, C20 ~ C20 ~ Clo + C20 + Clo ~ ^(20 + 1, "1), 

where -C20 + Clo + Clo - Clo " Clo + 2C20 - 1 = 0.2840... . 

For example. 



1 ^ 

1 

-1 
-1 

V 1 ) 




v 


+1 = - 







-2CI - 2CI 



2(5^ + + 1 

2C| + 2CI + 1 



(2C| + 2Cl + l)L 



\ -2Ci - 2Ci - 1 / 



In fact, 2(5^+2(5^ + 1 = -V5 ^ -2.236.... In other words, 
is an eigenvector of with eigenvalue X = —y/b. 

The SAGE code used to help with this calculation: 



(0,1,-1,-1,1) 



I SAGE Version 1.7.1, Release Date: 2007-01-18 I 
I Type notebookO for the GUI, and licenseC) for information. I 



sage : quadratic^residues (5) 
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[0, 1. 4] 

sage : quadrat ic_residues (7) 
[0, 1, 2, 4] 

sage: MS = MatrixSpaceCCyclotomicFieldCS) ,5,5) 
sage: V = VectorSpaceCCyclotomicFieldCB) ,5) 
sage: v = ¥([0,1,-1,-1,1]) 
sage: z = CyclotoniicField(5) ,genC) 

sage: F5 = MS( [[1, 1, 1, 1, 1] , [l,z,z-2,z-3,z-4] , [1 ,z-2,z-4,z~6 ,z~8] , [1 ,z~3 ,z~6,z-9,z- (12) ] , [1 ,z-4,z-8,z- (12) ,z- (16) ] ] ) 
sage: latex CF5) 

\lef t (\begin-[array>{rrrrr} 
l&l&l&l&lW 

Iazeta5&zeta5"{2}azeta5"{3}a-zeta5~{3} - zeta5"-[2} - zetaS - 1\\ 
l&zeta5"-[2}&-zeta5"{3} - zeta5"{2} - zeta5 - I&zeta5&zeta5"{3}\\ 
l&zeta5"-[3}&zeta5&-zeta5"-C3} - zeta5"{2} - zeta5 - I&zeta5"-C2}\\ 
l&-zeta5"{3> - zeta5"{2} - zeta5 - I&zeta5"{3}&zeta5'"{2}&zeta5 

\end{array}\right) 
sage : F5*v 

(0, -2*zeta5~3 - 2*zeta5'2 - 1, 2*zeta5"3 + 2*zeta5"2 + 1, 2*zeta5"3 + 2*zeta5"2 + 1, -2*zeta5"3 - 2*zeta5"2 - 1) 
sage: a = 2*z'3 + 2*z"2 + 1 
sage: a" 2 
5 

sage : exp(pi*I) 

-1.00000000000000 + 0.00000000000000323829504877970*1 

sage : zz=exp(2*pi*I/5) 

sage: aa = 2*zz"3 + 2*zz"2 + 1 

sage: aa 

-2.23606797749979 + 0.0000000000000581756864903582*1 

sage : zz=exp(-2*pi*I/5) 

sage: aa - 2*zz~3 + 2*zz"2 + 1 

sage : aa 

-2.23606797749979 - 0.0000000000000588418203051332*1 
sage: MS = MatrixSpace(CyclotoniicField(5) ,5,5) 
sage: z = CyclotomicField(5) -genC) 

sage: F5 = MS( [[1, 1, 1, 1, 1] , [l,z,z"2,z"3,z~4] , [l,z~2,z"4,z"6,z"8] , [1 ,z-3,z"6,z"9,z" (12)] , [1 ,z"4,z-8,z" (12) ,z" (16)]] ) 
sage: F5,f cp() 

(x + -2*zeta5~3 - 2*zeta5"2 - 1) * (x + 2*zeta5"3 + 2*zeta5"2 + 1)~2 * Cx"2 + 5) 
sage: F5.f cp() 

(x + -2*zeta5'3 - 2*zeta5~2 - 1) * (x + 2*zeta5~3 + 2*zeta5"2 + 1)"2 * (x"2 + 5) 
sage: z"4 

-zeta5"3 - zeta5"2 - zeta5 - 1 
sage : MS , identity_matrix () 

[1 0] 
[0 1 0] 
[0 1 0] 

[0 1 0] 

[0 1] 

sage: A = F5 - C-2*z~3 - 2*z'2 - 1) *MS . identity_matrix() 
sage : A . kernel ( ) 

Vector space of degree 5 and dimension 2 over Cyclotomic Field of order 5 and degree 4 
Basis matrix: 

[1 -zeta5"3 - zeta5''2 - 1 -zeta5"3 - zeta5"2 - 1 0] 

[0 1-1-1 1] 

sage: A.kemelC) .basisO 

[ 

(1, 0, -zetaS^S - zeta5"2 - 1, -zeta5"3 - zeta5"2 - 1, 0), 

(0, 1, -1, -1, 1) 
] 

sage: -z"3 - z"2 - l==z+z"4 
True 

sage: A = F5 + C-2*z"3 - 2*z"2 - l)*MS.identity_matrix() 
sage: A.kemelO .basisO 

[ 

(1, l/2*zeta5"3 + l/2*zeta5"2, l/2*zeta5"3 + l/2*zeta5"2, l/2*zeta5"3 + 

l/2*zeta5"2, l/2*zeta5"3 + l/2*zeta5'2) 

] 

sage: zz+zz"4 

0.618033988749874 + 0.0000000000000115463194561016*1 
sage : 4* (zz+zz"4) 

2.47213595499949 + 0.0000000000000461852778244065*1 
sage: Czz+zz"4) "2 

0.381966011250079 + 0.0000000000000142720357376695*1 
sage : (zz+zz'"4) "4 
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0.145898033750295 + 0.0000000000000109028651262724*1 
sage : (zz+zz''4) "5 

0.0901699437494591 + 0.00000000000000842292652849006*1 
sage : Czz+zz"4+l) "2 

2.61803398874982 + 0.0000000000000373646746498727*1 

sage: l/2*zz"3 + l/2*zz~2 

-0.809016994374949 - 0.0000000000000147104550762833*1 

sage: MSI = MatrixSpace(CyclotomicFieldC20) ,5,5) 

sage: z20 = CyclotomicFieldC20) . genC) 

sage: z5 = z20''4 

sage: I in CyclotomicFieldC20) 

False 

sage: z4 = z20"5 

sage: z4"4 

1 

sage : z4 
zeta20"5 
sage: z4"3 
-zeta20"5 

sage: F5I = MS ( [ [1 , 1 , 1 , 1 , 1] , [1 , z5 , z5"2 , z5"3 ,z5"4] , [1 ,z5"2,z5"4 ,z5"6 ,z5"8] , [I,z5'3,z5"6,z5"9,z5"(12)] , [l,z5"4,z5"8,z5'Cl2) ,z5"Cl6)]]) 
sage: A = F5I + z4*C-2*z5"3 - 2*z5~2 - l)*MSI.identity_matrixC) 
sage: A.kemelC) .basisC) 

[ 

CO, 1, zeta20"7 + zeta20"6 - zeta20"5 - zeta20"4 + zeta20"3 - 2*zeta20 - 1, 

-zeta20"7 - zeta20"6 + zeta20"5 + zeta20"4 - zeta20"3 + 2*zeta20 +1,-1) 

] 

sage: a = z20"7 + z20~6 - z20"5 - z20"4 + z20"3 - 2*z20 - 1 

sage: a" 5 

165*zeta20"7 + 125*zeta20"6 - 105*zeta20"5 - 125*zeta20"4 + 45*zeta20"3 - 210*zeta20 - 193 

sage: a in CyclotomicFieldC5) 

False 

sage: zz20=expC-2*pi*I/20) 
sage: zz20"5 

0.00000000000000127675647831893 - 1.00000000000000*1 

sage: A = F5I - z4*(-2*z5~3 - 2*z5"2 - l)*MSI.identity_matrixC) 

sage: A.kemelC) .basisC) 

[ 

Co, 1, -zeta20"7 + zeta20"6 + zeta20"5 - zeta20"4 - zeta20*3 + 2*zeta20 - 1, 

zetaSO"? - zeta20'6 - zetaSCB + zeta20~4 + zeta20"3 - 2*zeta20 + 1, -1) 

] 

sage: -zz"3 - zz"2 - 1 

0.618033988749899 + C.CCCCCCCCCCC00294209101525666»I 
sage: l/2*zz"3 + l/2*zz"2 

-0.809016994374949 - 0. 0000000000000147104B50762333H 

sage: -zz20"7 + zz20"6 + zz20"5 - zz20~4 - zz20"3 + 2*zz20 - 1 

0.284079043840412 + . 000000000000000222044604925031«I 

sage: zz20'7 + zz20'6 - zz20"5 - zz20"4 + zz20"3 - 2*zz20 - 1 

-3.B2014702134020 - 0.00000000000000222044604926031*1 

sage: [j for j in rangeCB)] 

[0, 1, 2, 3, 4] 

sage: [j*2y,5 for j in rangeCS)] 
[0, 2, 4, 1, 3] 

sage: [j*3°/,5 for j in rangeCS)] 
[0, 3, 1, 4, 2] 

sage: [j*4%5 for j in rangeCB)] 
[0, 4, 3, 2, i: 
sage; FB.fcpC) 

(x + -2*zeta5"3 - 2*zeta5"2 - 1) * Cx + 2*zeta5"3 + 2+zeta5"2 + 1)"2 * Cx"2 + 5) 

sage: -2*zz"3 - 2*zz"2 - 1 

2.23606797749979 + 0.0000000000000588418203051332*1 
sage: 



The table below lists the multiplicity of the eigenvalues of F^, for some 
small values of N: 



58 





mult, of 


mult, of 


mult, of 


mult, of 


N 


^/n 




i^/N 


-i^/N 


4 


2 


1 





1 


5 


2 


1 


1 


1 


6 


2 


2 


1 


1 


7 


2 


2 


1 


2 


8 


3 


2 


1 


2 


9 


3 


2 


2 


2 


10 


3 


3 


2 


2 


11 


3 


3 


2 


3 


12 


4 


3 


2 


3 


13 


4 


3 


3 


3 



In general, the multiplicity of A = is equal to the rank of M^^^ in 

fjTijl . According to Good [HI, this is 

~m + 4)] [UN + 2)] m-l)] [iiN + W 
Here is the SAGE code verifying the last line of the first table above: 

sage : p = 13 

sage: MS = MatrixSpaceCCyclotomicFieldCp) ,p,p) 
sage: z = CyclotomicFieldCp) ,genC) 
sage: zz = exp C-2*pi*I/p) 

sage: r = lambda k: [z"Cj*k) for j in rangeCp)] 
sage: F = MSCErCk) for k in rangeCp)]) 
sage: F.fcpC) 

(x + -2*zetal3"ll - 2*zetal3"8 - 2*zetal3"7 - 2*zetal3"6 - 2*zetal3"5 - 2*zetal3"2 - 1)"3 * 

Cx + 2*zetal3"ll + 2*zetal3"8 + 2*zetal3~7 + 2*zetal3"6 + 2*zetal3"5 + 2*zetal3"2 + 1)"4 * Cx"2 + 13)"3 
sage: 2*zz"ll + 2*zz"8 + 2*zz"7 + 2*zz"6 + 2*zz"5 + 2*zz"2 + 1 
-3.60555127546399 - 0.00000000000000177635683940025*1 

Here is an example of such a function. Let p be a prime number and 
let £{j) = 1 if j is a non-zero square mod p and —1 is j is a non-zero non- 
square mod p and £{0) = 0. This is called the Legendre function or the 
quadratic residue symbol. It turns out that (since the product of two squares 
is a square and the product of a square with a non-square is a non-square) 
that i{jk) = £{j)£{k)^ for all non-zero j, k. Also, it can be checked that, for 
any < k < p^ the set of multiples of k mod p is the same as the set of all 
integers mod p: 

{jk I J e z/pZ}. 
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Therefore, the assumptions made in and above liold. The eigenval- 
ues are sometimes called Gauss sums but we shall not discuss them further. 
Further information are available in the excellent papers by Good and 
McClellan and Park pPj . 

3.2 The DFT and the coefficients of the FS 

We saw in fjlO|) the approximation which motivated the definition of the DFT. 
If / is a "nice" function on (0, P) and if / = (/(|P), /(ip), /(i^P)) 
is the vector of sampled values of / then 

FS{f), ^ j^DFTif),, (17) 

where k ranges over {0, 1, — 1}. Based on the approximation in (|1U|). 
one expects that the estimate (fTTjl is only good when k/N is "small". 

Example 25 Let f{x) = for < x < 1, extended periodically to M with 
period P = 1. The graph looks something like: 



Ck 



Figure 11: Graph of /(x), — 1 < x < 2. 
We compute its k-th Fourier series coefficient: 

io Jo Jo 



1 + 27ri/c' 



60 



for all k E Zi. Now, let 



f= (/o,/i,...,/7v-i) = (l,e^^/^, 



-2/N 



^{N-l)/N^ 



be the vector of sampled values. We compute its k-th DFT component: 



DFT{f), 



sr^N-l r -2-nikj/N 
Yf^~^ f,-j/N^-2nikj/N 



l-e" 



l_g(-l-2TTik)/N 



1-e- 



l_g(-l-2TTik)/N ■ 

The estimate (|i7| ) simply asserts in this case that 

1 1 - 1 - 



1 + 2mk 

Is this true? 

If we use the approximation 



I _ g(-l-27rifc)/Af ' 



e ^ = 1 - X + -X + ... 

we see that, for "small" (-1 - 2mk)/N, 1 - e(-i-2^*'=)/^ ^ ^(1 + 2mk). 
With N = 10, even the first few values aren't very close. 



k 


Ck 


DFTif), 


\ck-DFT{f)k\ 





0.63212 


0.664253 


0.03213 


1 


0.01561 - 0.0981 U 


0.04775 - 0.094781 


0.03231 


2 


0.003977 - 0.049981 


0.03615 - 0.043191 


0.03288 


3 


0.001774 - 0.03344i 


0.03401 - 0.022871 


0.03392 



Here is the list plot of the values of \ck — DFT{f)k\ 
and here is the histogram plot: 
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Figure 12: 


List plot of 


\ck-DFT{f)k\,k^0,...,9. 





You can see from these graphs that the larger k/N gets, the worse the 
approximation is. 

When N = 100 the approximation is about 10 times better, as is to be 
expected. 



k 


Ck 


DFT{f)k 


Cfc - DFT{fW 





0.6321 


0.6352 


0.003173 


1 


0.01561 - 0.09812*1 


0.01878 - 0.09808*1 


0.003165 


2 


0.003977 - 0.04998*1 


0.007143 - 0.04992*1 


0.003166 


3 


0.001774 - 0.03344% 


0.004939 - 0.03334*1 


0.003167 



Here is the list plot of the values oflO\ck — DFT{f)k\ (the error has been 
scaled up by a factor of 10 so that the plot comes out nicer): 
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0.3 
0.2 



-0.2 - 
-0.3 ; 



Figure 13: Histogram of \ck — DFT{f)k\^ /c = 0, ...,9. 



Figure 14: List plot of \ck - DFT{f)k\, k = 0,...,99. 
Here is the SAGE code for the = 10 case: 



sage: CC5 = ComplexField(15) 
sage : N = 10 
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sage: ck - lEonbda k; Cl-expC-1) )/Cl+2*pi*I*k) 

sage: dftk - lambda k: (l/N)*(l-exp(-l))/(l-e3:pC (-l-2*pl*I»k)/N)) 
sage: [CC5Cck(J)) for j in rangsCN)] 

[0.6321, 

0.01561 - 0.09812*1, 
0.003977 - 0.04998*1, 
0.001774 - 0.03344*1, 
0.0009991 - 0.02611*1, 
0.0006397 - 0.02010*1, 
0.0004444 - 0.01675*1, 
0.0003265 - 0.01436*1, 
0.0002500 - 0.01257*1, 
0.0001976 - 0.01117*1] 
sage; [CCBCdftkCj)) for j in rangeCN)] 

[0.6642, 

0.04775 - 0.09479*1, 
0.03615 - 0.04319*1, 
0.03401 - 0.02287*1, 
0.03334 - 0.01024*1, 
0.03318 - 0.00000000000000006104*1, 
0.03334 + 0.01024*1, 
0.03401 + 0.02287*1, 
0.03615 + 0.04318*1, 
0.04775 + 0.09478*1] 
sage: [abB(CC6(ck(j))-CC6(dftk(j))) for j in range(N)] 

[0.03213, 
0.03231, 
0.032SS, 
0.03392, 
0.03660, 
0.03825, 
0.04256, 
0.05021, 
0.06631, 
0.1161] 



sage: 


L " [abs(CC5(ckCj))-CC5(dftk(j))) for j in range(N)] 




sage: 


show ( list _plot(L) ) 




sage: 


J = range CN) 




sage: 


s = IndexedSequence (L , J) 




sage: 


(s .plot_histograinC)) . save C "hist ogram-dftk-vs-cklO.png" ,3cmin=-l,3ai 


iax=10,yinin=-0 .6 ,yma3:=0 .6) 



Let's try one more example. 

Example 26 Let f{x) — for < x < 1, extended periodically to M. with 

period P = 1. The graph looks something like: 

You can enter this function and plot it's values, for —l<x<2, in SAGE 
using the commands 

sage: fO = (x+l)-2; fl = x"2; f2 = (x-l)-2 

sage: f = Piecewise ( [ [(-1,0) ,f 0] , [(0, 1) ,f 1] , [(1 , 2) ,f 2] ] ) 

sage: show (f .plot () ) 

We compute, for k ^ 0, 



Here is the plot of the real part of the partial sum ^_io c^e' 
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■As ' 



Figure 15: Graph of f{x), —l<x<2. 
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Figure 16: Graph of Sio{x), -1 < x <2. 



sage: c = lambda k : I* ( C2*k"2*pi"2 - l)/(4*k"3*pi"3) + 1/ C4*k"3*pi"3) ) + 1/ C2*k"2*pi"2) 
sage: ps = lambda x: 1/3+sumC [CcCk)*expC2*I*pi*k*x) ) .realC) for k in rangeC-10,10) if k!=0]) 
sage : showCplot (ps ,-1,2)) 
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3.3 The DFT and convolutions 



This section is based, not on Walker's book |Wlj but on material in Frazier's 
well- written book [Fj. 

Let 'L/N'L denote the abelian group of integers mod A^, and let 

VN = {f\f: Z/iVZ ^ C}. 

This is a C-vector space which we can identify with the vector space via 
the map / i — > {f{0), /(I), f{N-l)), Vn ^ C^. You can visualize Z/NZ 
as a circle with N equally spaced points and functions on Z/NZ as "weights" 
on each of these points: 
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Figure 17: Circle representing Z/iVZ. 
Define convolution by 

if, 9) ^ f*9, ^ ^ 

where 

(f*gm^ J2 m9(k-e). 

— * — * 

In other words, if / = (/o, /i, /iv-i) and g = (g^o, 9i, 9n-i) then / * ^ is 
another vector, whose k-th coordinate is given by 

(/ * 9)k = f^9k-e, 

where are subscripts are computed mod and represented in the set {0, 1, A^— 
1} This binary operation on Vn is commutative. In other words, f *9 — 9* f' 
ife'^k-e then 
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We will be done once we show that (for each k e Z/A^Z) as ^ ranges over 
all of 'L/N'L, so does £' — k — i. But if i' misses something in Z/NZ, say x, 
then k — £ X ioT all £ E Z/NZ, and so £ ^ k — x. This is a contradiction, 
so Etez/ATz f{^)9{k -£)^ ileez/Nz /(^ - -^Of (^') = * as desired. 

Example 27 Here is a SAGE session for computing the convolution: 

sage: F = CyclotoiiiicField(16) 
sage: J = range (16) 

sage: A = [F(0) for i in J] ; A[0] = F(l); A[l] = F(l); A[2] = F(l); A[3] = F(l) 

sage: z = F.genO 

sage: B = [z~(i) for i in J] 

sage: sA = IndexedSequence(A, J) 

sage: sB = IndexedSequence(B, J) 

sage: sA. convolution(sB) 

[1. 



zetaie + 1, 








zetal6"2 + 


zetaie + 1, 






zetal6"3 + 


zetal6"2 + 


zetaie + 1, 




zetal6''4 + 


zetal6~3 + 


zetaie"2 + 


zetaie , 


zetaie'S + 


zetal6~4 + 


zetaie"3 + 


zetaie"2. 


zetal6"6 + 


zetal6"5 + 


zetal6"4 + 


zetal6"3, 


zetal6"7 + 


zetal6"6 + 


zetaie"5 + 


zetaie"4, 


zetaie"? + 


zetal6~6 + 


zetaie"5 - 


1, 


zetal6~7 + 


zetal6~6 - 


zetaie - 1, 




zetal6"7 - 


zetal6"2 - 


zetaie - 1, 




-zetal6"3 - 


• zetaie"2 - 


■ zetaie - 1, 


-zetal6''4 - 


• zetaie"3 - 


zetaie"2 - 


zetaie. 


-zetaie^S - 


• zetaie"4 - 


zetaie~3 - 


zetaie"2, 


-zetal6"6 - 


■ zetal6"5 - 


• zetaie"4 - 


zetal6"3, 


-zetal6"7 - 


• zetaie"e - 


■ zetaie"5 - 


zetal6"4, 


-zetaie"? - 


• zetaie"e - 


zetaie"5. 




-zetaie"7 - 


• zetaie"e. 






-zetal6"7, 








0, 0, 0, 0, 


0, 0, 0, 0, 0, 0, 0, 


0] 
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Let ( = denote a primitive A^*^ root of unity in F (W, in the notation 
above, is one good choice). Recall, for g e Vn, the discrete Fourier transform 
of g was defined by 

^^(A)= Yl AeZ/TVZ. 
Also, we defined the inverse discrete Fourier transform of G by 

A basic and very useful fact about the Fourier transform is that the 
Fourier transform of a convolution is the product of the Fourier transforms. 
Here's the proof: 



(/ * = Eeem^ ^>^ezm /(^)^(^ - 

1,/NZ /( 

r(A)5^(A) 



In coordinate notation: 

DFT{f* g)k = DFT{f\DFT{g)k, 
ior al\0<k<N - 1. 

Definition 28 For /i G V^jv, define Mh : Vn ^ Vn by 

If X denotes the componentwise product of two vectors, 

(oo, ai, a^v-i) X (60, ^at-i) = {aobo,aibi, ...,aN-ibN-i), 
then, in coordinate notation, 

Mj;if)^DFT-\hxDFTif)), 

for f,hE ■ A linear transformation T : Vat ^ Vat of the form T = M/j, 
for some h E Vn, is called a Fourier multiplier operator. 



69 



In other words, a Fourier multiplier operator (represented in the standard 
basis) is a linear transformation of the form F^^DFn, where D is an N x N 
diagonal matrix. Note that the product of two Fourier multiplier operators 
is a Fourier multiplier operator: Mh^M^,^ = M^^h.^. 

Stereo systems and FMOs: Here is one way in which Fourier multiplier oper- 
ators can be thought of in terms of Dolby stereo (this is a grossly over-simplied 
description but you will get the idea). Dolby cuts off the high frequencies, often 
which are crackles and pops and other noise in the channel, making the music 
sound nicer. How does it do that? If / G Vn represents the digital sample of the 
sound, DFT{f) = represents the frequencies of the sound. To cut off the high 
frequencies, multiply DFT{f) by some ("low-pass filter") h G Vn which is on 
the high frequencies and 1 on the rest: you get hf^. To recover the sound from 
this, take the inverse DFT, (hf^)"^. This is the same sound, but without the high 
frequencies. This "filtered sound" is an example of a Fourier multiplier operator. 

The following result appears as Theorem 2.19 of jFj. It characterizes the 
Fourier multiplier operators. 

Theorem 29 Let T : Vn Vn denote a linear operator. The following are 
equivalent: 



1. 


T is 


translation invariant. 




2. 


The 


matrix A representing T in the 


standard basis is circulant. 


3. 


T is 


a convolution operator. 




I 


T is 


a Fourier multiplier operator. 




5. 


The 


matrix B representing T in the 


Fourier basis is diagonal. 



We shall define all these terms (convolution operator, etc) give some ex- 
amples, and prove this theorem. 

Remark 1 As a consequence of the proof below, we shall show that, for all f,g€ 
Vn, 

Let M{N) denote the number of multiplications required to compute the DFT on 
Vjsf. The above identity implies that the number of multiplications required to 
compute the convolution / * is at most 2 • M{N) + 1. We shall see in ^6.11 that 
M{N) < iV(log2(iV) + 2). (This remark shah be applied in below.) 
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Definition 30 • Let r : V^v Vn denote the translation map: {Tf){x) = 
f{x + 1) (addition in Z/NZ). Note 



r : (xo 



{Xi,X2, ...,Xn-1,Xo). 



• We say an operator T is translation invariant if the diagram 



T 



N 



T 



T 



Vn ^ Vn 

commutes. (This phrase will be explained below.) 

• Define the convolution operato^ associated to g, 



hjT,{f) = f*g. 



Remark 2 Here is a remark on the grammar used in the diagramatical 
definition of translation invariance above. The phrase "diagram commutes" 
is a fancy way to say that, for each / G Vat (picking an element in the copy 
of Vjy in the upper left hand corner), the element T(r(/)) G Vat (mapping 
from the upper left corner along the top arrow and down the right arrow 
r(/) I — >• T(r(/))) is equal to the element r(T(/)) (mapping down the left 
arrow and along the bottom arrow), as functions on Z/A^Z. In other words, 
T is translation invariant if and only if, for all k G Z/NZ and all f E Vn, 



Example 31 Let's look again at the operator neg : Vn — > Vn, which sends 
f{k) to f{—k). Is this translation invariant? To answer this, we must see 
whether or not (r(neg/))(A;) = (neg(r/))(A;), for each k G Z/NZ and each 



/ G Vn- We have, for example, (r(neg/))(0) = (neg/))(l) = /(—I) = 
f{N-l) and (neg(r/))(0) = (r/)(-0) = (r/)(0) = /(I). In general, these 



two coordinares are different, so neg is not translation invarianty. 

^As usual, the term "operator" is reserved for a linear transformation from a vector 
space to itself. 

^ However, if you restrict neg to the subspace of Vn of functions / for which f{N — £) = 
f{£) then it is translation invariant (in fact, on this subspace, it is the identity). 



we have T{T{f)){k) 
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Lemma 32 The convolution operator Tg is translation invariant. In other 
words, the diagram 



Vj 



N 



N 



V 



N 



1". 



commutes, for all g e Vn- 

In terms of the above theorem, this lemma says "3 
proof: Recall 

for k e Z/NZ. We have 



Wf)m ={T{f)*g){k) 

= j:iem^irfmg{k-i) 

^Tg{f){k + l)^T{Tg{f)){k). 



□ 



Definition 33 An N x N matrix A is circulant if and only if 

^k,e — ^fc+l(mod N), £+l(mod N), 

for aU < < iV - 1, < £ < iV - 1. 

In other words, to go from one row to the next in a circulant matrix, just 
"rotate" or cyclically permute the rows. In particular, (setting k = £ = 1) 
Al l = ^2,2 = ••• = ^N,N, SO all the diagonal entries must be the same. 

Example 34 Let N — 5. The matrix 
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/ 1 

5 
4 
3 
V2 



3 4 5 \ 
2 3 4 
1 2 3 
5 1 2 

4 5 1/ 



is circulant. 



Lemma 35 A linear transformation T : Vn Vn is translation invariant 
if and only if the matrix representing it in the standard basis is circulant. 

In terms of the above theorem, this lemma says "1 2" . 

proof: Since T : Vjy ^ Vj^ is hnear, with respect to the standard vasis it 
is represented by an x A*" matrix 

Tf^Af, A={Aj), 

where / = *(/(0), /(I), /(iV-1)). In other words, T(/)(fc) = Eeez/Nz^k/fii), 
for k e Z/NZ. For k e Z/NZ, we have 

T{T{f)){k) =Eee^/^zAkArfm 
= 'E.iez/NZ ^k,e.f{( + 1) 
= Y.e'ez/NZ ^k/'~if{^'), 

where the 2""' subscript of A^j is taken mod A^. Changing the dummy variable 
£' to i, this becomes 

T{r{f)){k)= J2 ^M-i/W- 

On the other hand, 

riT{f)){k)= J2 

Comparing coefficients, it follows that the linear transformation T is transla- 
tion invariant if and only if its matrix A satisfies: A^^^ — A^^i (mod n), i+i (mod n), 
for aR < k < N - 1, < e < N - 1. But this is true if and only if A is 
circulant. □ 
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Example 36 Let = 5 and let T : ^ be defined by 

T{f){k) = {CfX 

where C is as in Example EH and / is the column vector / = (/(O), /(I), /(2), /(3), /(4)), 
for / G V5. Here is an example: 



k 


12 3 4 


m 


10 10 


rf{k) 


10 10 


Tf{k) 


4 7 5 8 6 




6 4 7 5 8 


nrfm 


6 4 7 5 8 



We see in this example that Tr^f) = TT{f), consistent with the fact 
"(2) =^ (1)". 

Lemma 37 // a linear transformation T : is translation invariant 

then it is a convolution map. 

In terms of the above theorem, this lemma says "1 =^ 3" . 

proof: Let T denote a translation invariant linear operator. Define g G 
Vn by g{k) = Ao_k (mod n), k e Z/NZ. Then 

nfm= Yl ^o,if{i)= E 9{-i)m, 

for all / G Vat. Replacing / by a translation (mod A^) (note g is periodic 
with period A^) gives 

= j:e'ez/Nz9{k-nf{i') {£' = i + k) 
= J2eez/Nz9ik — i)fii), 

for all / G Vat. In other words, T is a convolution map. □ 

Using the above lemmas, we see the connection between maps given by 
circulant matrices and convolution operators. 
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Lemma 38 IfT: Vn —>■ Vn is a convolution operator then, represented in 
the Fourier basis, T is diagonal. 

This says (3) =^ (5). 

proof: By hypothesis, T = Tg, for some g E Vn- Let us compute the 
action of Tg on the function 6^ G Vn, where bk{i) = e^'^*'^^/^. We have 

T,mk) =j:e^^/N^e'-'^'^^gik-i) 

so Tg(hj) = g^bj. If T is represented hj A = {Aj^k), in the Fourier basis, then 

T{bj) = ^Aj^kh, 
j 

for < J < iV - 1. Since the {bj}rQ are hnearly independent (they are a 
basis, after all!), we may compare coefficients to conclude: Ajj = g^{j) for 
all j and Aj^k = ii j ^ k. In other words, A is diagonal, as desired. □ 
In terms of the above theorem, this computation proves "3 =^ 5" . 

Lemma 39 T : V^v — > V^v is a convolution operator if and only if T is a 
Fourier multiplier operator. 

This says (3) (4). 

proof: Suppose T : V^v — Vat is a convolution operator, say T = Tg, 
for some g G Vn- Let h = g^ denote the inverse discrete Fourier transform 
of g. The claim is that Tg = M^, where is defined as in Definition 1281 
Since the Fourier transform of the convolution is the product of the Fourier 
transforms, for each / G Vn, we have {f *g)^ = f^g^- Taking inverse Fourier 
transforms of both sides gives 

Tg{f) = f*g=ig^fr=Mg.{f). 

Therefore, if T is a convolution operator then it is also a Fourier multiplier 
operator. 

Conversely, suppose T is a Fourier multiplier operator, say T = M^, for 
some h G Vn- Let g = - The claim is that Tg = Mh- The proof of this 
case also follows from the identity (/ * g)'^ = f^g^- □ 
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Lemma 40 If T : Vjv —>■ Vn is diagonal, represented in the Fourier basis, 
then T is translation invariant. 

This says (5) =^ (1). 

proof: By hypothesis, there are constants 7^ G C such that T{be){k) = 
^(bi{k), for each C. with < £ < — 1 and each k G Z/A^Z. Therefore, if 
/ G Vat is expressed in terms of the Fourier basis as f\x) = (for 
Ci G C) then 

T{f){k) = Y^caMk). k G Z/iVZ. 

e 

Now we compute 



r(T/)(A;) = T/(A; + 1) = Y^^aMk + 1) = e-2-^/^Q7A(fc), 

and 

r/(fc) = Y.c,Th,{k) = Y^cMk + l) = Y,e-^-'"''cMk). 
i I I 

so 

£ £ 

This imphes TT{f) = Tr^f), for all / G Vat, as desired. □ 

We have: (3) =^ (1) (LemmaEl, (1) ^ (2) (LemmaESl), (1) ^ 
(3) (Lemma E7|), (3) =^ (5) (Lemma |^, (3) (4) (Lemma EHD, and 

(5) =^ (1) (Lemma l4Up. These lemmas taken together finished the proof 
of the theorem. 

In the next section, we shall give an example of how the Cesaro filter 
gives rise to a Fourier multiplier operator. 

4 Filters and reconstruction 

Here is a problem which is sometimes called the reconstruction problem: Sup- 
pose we know the Fourier series coefficients c„ of f{x) but we don't know 
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f{x) itself. How do we "reconstruct" f{x) from its FS coefficients? The the- 
oretical answer to this is very easy, it is (by Dirichlet's Theorem ^ simply 
the FS expansion: 



Suppose that we modify this problem into a more practical one: Suppose 
we know the Fourier series coefficients c„ of f{x) but we don't know f{x) 
itself. Given some error tolerance 6 > 0, how do we "reconstruct" (efficiently) 
f{x), with error at most 6, from its FS coefficients? Practically speaking, a 
"filter" is a sequence of weights you apply to the FS coefficients to accomplish 
some aim (such as approximating f{x) "quickly", or filtering out "noise", or 
...). Can we weight the terms in the partial sums of the FS to compute 
an approximation to the value of a FS in a more efficient way than simply 
looking at partial sums alone? In this section, we examine several filters and 
see that the answer to this question is, in a specific sense, "yes" . 

4.1 Dirichlet's kernel 

Let 



denote the M-th partial sum of the FS of /. (This is a filter, where the 
weights are all 1 from —M to M and elsewhere. Sometimes the term 
"rectangular window" is seen in the literature for this.) Here is an integral 
representation for Sm- First, recall that 



M 



SMix) = 



,2iTijx/P 



j=-M 




SO 



Sm{x) 



J,''f{t)K^ix-t)dt 
f*Kg{x), 



(19) 
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where K^{z) = ^J^jL-m^'^'^^'^^^^- (Note that Walker defined the convo- 
lution in a slightly different way than we do.) This function is called the 
Dirichlet kernel . It turns out this function has a nice closed-form expres- 
sion: 

Some plots of this, for various values of M, are given below. You can see that 
the graphs get 'spikier and spikier" (approaching the Dirac delta function, 
5) as M gets larger and larger. 




sage: M = 5 

sage: f = lambda z: (l/(M+l))*(sin((2*M+l)*z/2)/siii(z/2)) 

sage: PI = plot (f ,-5,5) 

sage: M = 10 

sage: f = lambda z: (l/(M+l))*(sin((2*M+l)*z/2)/sin(z/2)) 

sage: P2 = plot (f ,-5,5) 

sage : M = 50 

sage: f = lambda z: (l/(M+l))*(sin((2*M+l)*z/2)/sin(z/2)) 
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sage: P3 = plot (f ,-5,5) 
sage: show(Pl+P2+P3) 



Here is the proof of ()2()j) : 



27Tijz/P 



as desired. 



„-2mMz/P sr^2M 2nijz/P 
^-2-KiMzlP e'^^^^^M+l)z/P_^ 

^ ^2iTiz/ P _l 

^-2niMz/P {e^'<-^^'+''>^/P_ 



-7ri(2M+l)z/P 



)/2i 



(e7riz/P_e-Tiz/p-)/2j 
^-2niMz/P+TTi(2M+l]z/P~TTiz/P sin{n{2M+l)z/P) 

sin{iTz/P) 

sin(7r(2M+l)z/P) 
sin(7r^/P) ' 



4.2 Cesaro filters 

Aside: Here is a historical/mathematical remark explaining the idea behind this. Start 
with any infinite series, 



do + fll + ^2 + 03 + •• 



which, for the sake of this discussion, let's assume converges absolutely. This means that 
the series of partial sums 



So = flo, si = flo + ai, , S2 = flo + ai + 02, ... 

has a limit - namely the value of the series aj. In particular, the series of arithmetic 
means 



m2 



Too = So = 
3 



flo, mi = = i(ao + ao + fli) = flo + _ 

i(ao + flo + ai + flo + fli + 02) = flo + lai + 



2^1' 



-.0.2, 



also has a limiting value, namely the value of the series aj. (After all, you are simply 
averaging values which themselves have a limiting value.) In general. 



TOJ 



J 

E 

i=i 



(1- 



J+1 



In the early 1900's the basic observation that the limit of the mj's, as J ^ 00, is equal to 
ttj was applied by Fejer to the study of Fourier series. A more general convergence test 
was devised by Cesaro, who generalized the "weights" 1 — used above. For historical 
reasons, the "weights" 1 — are sometimes called Cesaro filters. □ 
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As usual, let Sm{x) = ^Y^j=-M^^^^^^^^^ denote the M-th partial sum of 
the FS of / and let 



The above historical discussion motivates the following terminology. 

Definition 41 The M-th Cesdro-filtered partial sum of the FS of / is the 
function 5*^^ = Sf.^j above. 

As the graphs show, as M gets larger, the Sf/s approximate / "more 
smoothly" than the S'a/'s do. The factor (1 — -^) have the effect of "smoothing 
out" the "Gibbs phenomenon spikes" you see near the jump discontinuities. 

Example 42 As an example, here are the plots of some partial sums of the 
Fourier series, and filtered partial sums of the Fourier series. 
Let f{x) = e""*, < X < 1 as in Example above. 

We shall use list plots, since they are easy to construct in SAGE . Here is 
f again, hut as a list plot: 




j=-M 



j=0 



Figure 19: List plot of values ofe ^,0<x<l. 
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Sr. 



Figure 20: List plot of values of £'5(2;), < x < 1. 



Figure 21: List plot of values oi {x), < x < 1. 



Note how the Gibbs phenomenom of is "smoothed out" by the filter 
the graph of seems to be less "bumpy". 
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Figure 22: List plot of values of S'io(a;), < x < 1. 



Figure 23: List plot of values of Siq{x), < a; < 1. 
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25- 



Figure 24: List plot of values of S25{x), < x < 1. 



^25' 



Figure 25: List plot of values of 5'£(a;), < a; < 1. 



In each case, we see that the Cesaro filter "smooths out" the the Gibbs 
phenomenom of the partial sums of the Fourier series. 
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Here is the SAGE code for the M — 25 case: 



sage: FSf = lambda x: CsuiiiC[ckCj)*expC2*pi*I*j*x) for j in raiigeC-25,25)] ) ) . realC) 

sage: L = [FSf Cj/50) for j in rangeCBO)] 

sage : showClist_plot (L) ) 

sage: FSf = lambda x: (sum( [Cl-absCj)/25)*ckCj)*exp(2*pi*I*j*x) for j in rangeC-25,25)] )) .realO 

sage: L = [FSf Cj/50) for j in rangeCSO)] 

sage : showClist_plot (L) ) 



Here is an integral representation for S^. First, recall that 

^ £ me-'-'^'/^ dt, 



so 



= ^ lo fit) E^-m(1 - g)e^-^(-*)/^ dt (21) 



= f{t)KM{x-t) dt, 



where 



P ^ ^ M' 

j=-M 

This function is called the Fejer kernel (it is also sometimes referred to as 
the "point spread function" of the filter). This has a simpler expression, 

1 1 ^ sin((M + l)7r2:/P) ^, 

This is included here not because we need it as much as because this expres- 
sion is much easier to graph: 

You sec how these functions seem to be, as M oo, approaching the 
spiky-looking Dirac delta function. In fact, (as distributions) they do. (A dis- 
tribution is a linear functional on the vector space of all compactly supported 
infinitely differentiable functions C^(M).) In other words, 

lim [ f{x)KM{x-t)dx^ I f{x)5{x-t)dx^f{t). 
Jr Jr 

This is the essential reason why, if / is continuous at x, limM^oo 'S'm(^) ~ 

/(^)- 
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Figure 26: List plot of values of 2ttKm{x), M — 5, 10, 50, P — 2tt (normalized 
by 27r for simplicity). 

The SAGE command for this: 

sage: M = 5 

sage: f = lambda z: (l/(M+l))*(sin((M+l)*z/2)/sin(z/2))-2 

sage: PI = plot (f ,-5,5) 

sage: M = 10 

sage: f = lambda z: (l/(M+l))*(sin((M+l)*z/2)/sin(z/2))-2 

sage: P2 = plot (f, -5, 5) 

sage: M = 50 

sage: f = lambda z: (l/(M+l))*(sin((M+l)*z/2)/sin(z/2))-2 

sage: P3 = plot (f, -5, 5) 

sage: show(Pl+P2+P3) 

4.3 The discrete analog 

There is a discrete analog of the Cm — Cmj which may be regarded as a 
Fourier multiplier operator. This subsection is included for the purpose of il- 
lustrating the notion of a Fourier multiplier operator by means of an example 
(as you will see, our process of discretizing the filter ruins its usefulness). 
If we replace the usual FT on (0, P) 
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/ 




by the DFT on C^, 



/ ^ DFT{f)k 



then the Cesaro filter becomes simply the modification 



k 
N 



)DFT{f\. 



If we define ^ G Vat by C,{k) = {1 — j^) then, in functional notation, the 
Cesaro filter becomes the modification 



This operator C is, by construction, a Fourier multiplier operator: C = M^. 
Here is an example computation. 

Example 43 Let f{x) = 10x(l — x), < a; < 1, which we sample at N = 25 
regularly spaced points. The plot of this function is an inverted parabola, 
passing through and 1 on the x-axis. Below, we plot both f and the real 
part ofC{f): 



Figure 27: List plot of values of 10x(l — x) and (the real part of) its image 
under a FMO. 




(/ e Vn). 

Vn by 
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4.4 Hann filter 

As usual, let S'm(x) = ^Y^=-M^k^^''^^^^ denote the M-th partial sum of the 
FS of / and let 

M 
j=-M 

where HMix) = (1 + cos(ff ))/2 is the Hann filte^ 

Definition 44 The M-th Hann-filtered partial sum of the FS of / is the 
function S^^ = S^^j above. 

Here is an integral representation for 5*1^. Since 




we have 

SUx) = E-L-M HMimio /(t)e-^-^^*/^ dt)e^-^^^IP 

= ^ lo m Ef=-M HM{j)e'-'^'^^-'y^ dt (22) 
= J^''fit)Kf,ix-t)dt, 

where 

M 
j=~M 

This function is called the Hann kernel. This has a simpler expression, 
where K?, is the Dirichlet kernel. 



®With 0.54 + 0.46cos(^) instead of i?M, you get the Hamming filter. We shall not 
describe this due to its similarity with the Hann filter. 
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Example 45 Consider the odd function of period 27r defined by 

f(^^ _ / -1> < X < 7r/2, 
' [2, tt/2<x< it. 

We use SAGE to compare the ordinary partial sum 520(2;) to the Hann- filtered 
partial sum S2o{x) and the Cesaro-filtered partial sum S^. As you can see from 
the graphs, the filtered partial sum is "smoother" than 5*20 (2;) and a slightly better 
fit than S2Q- 






Figure 28: Plot of f{x), S§^{x). Plot of f{x), Sg{x) 



Plot of /(x), S2o{x). 
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The SAGE command for this: 

sage: fl = lambda x:-2 

sage: f2 = lambda x:l 

sage: f3 = lambda x:-l 

sage : f 4 = lambda x : 2 

sage: f = Piecewise ( [ [(-pi , -pi/2) , f 1] , [(-pi/2 ,0) , f 2] , [(0 ,pi/2) ,f 3] , [(pi/2,pi) ,f 4] ] ) 

sage: PI = f .plot_f ourier_series_paxtial_sum_hEinn(20,pi , -5 ,5) 

sage: P2 = f .plot (rgbcolor= (0.8, 0.3, 0.4)) 

sage: P3 = f .plot_f ourier_series_partial_sum(20,pi,-5,5) 

sage: P4 = f .plot_fourier_series_pEirtial_sum_cesaro(20,pi,-5,5) 

sage: show(Pl+P2+P3+P4) 

4.5 Poisson summation formula 

Let S{x) be the Dirac delta function. The Poisson summation formula can 
be regarded as formula for the FT of the "Dirac comb" : 



We shall not need this formulation, or any of the many very interesting 
generalizations of this formula. In the next section, it will be applied to 
proving the Shannon sampling theory, so we only present here what we need 
for that. 

Theorem 46 (Poisson summation) Assume that P > is fixed and 

• f be a continuous function, 

• the function 



oo 



A{x) = S{x-Pn). 



n=—oo 



/p(^) = E/(^-p) 



converges uniformly on [—P/2, P/2\, and 




converges. 
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Then 



proof: This theorem says that the n — th Fourier series coefficient of the 
perioic function fp is 

Cn — pj\ p), 

where here / denotes the Fourier transform (on M). This identity is verified 
by the following computation: 



as desired. □ 



4.6 Shannon's sampling theorem 

Let / be a continuous function periodic with period P. So far, we have 
started with sampling values of f{x) at equally spaces points to get a 
vector / = (/(;^P))j=o,i,...,iv-i, then computed its DFT, which we showed 
was a good approximation to the FS coefficients: 

f\x) ^ f 

Ck = j; Jo /(x)e-2-^-/^ dx « DFT{f)k 

Based on this, we might expect that f{x) can be approximately reconstructed 

? 

from its sample values alone. For example, perhaps something like f{x) ~ 
^|j;.|^jy DFT(/)fee^'^*^^/^ might be true. Shannon's Sampling Theorem says 
that, under certain condtions, f{x) G L^(R) is completely determined from 
its sampled values. 
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We say / G L-^(M) is band limited if there is a number L > such that 
fit) = for all t with |t| > L. When such an L exists and is choosen as 
small as possible, the number 2L is called the Nyquist rate and the number 
2^ is the sampling period. 

Define the "sink" function 



sinc(x) — 




otherwise. 



Theorem 47 (Shannon's Sampling Theorem) Assume f is as in the above 
theorem and that f e L^(M) is band limited with Nyquist rate P = 2L. Then 

/(^) = /(^)si^c(2La; - n). 

neZ 



proof: Let fp be defined analogously to fp above. The Poisson formula 
gives 

My) -Enezhy-T) 

p L^n&L J \ pl^ 

— IV f(—ii\p'^'^iy/P 

since fix) = f{—x). By hypothesis, f{y) = for \y\ > P/2, so fp{y) is 
merely a periodic extension of /. Let 

= { 0, M > p/l: 

so, by hypothesis, f{y) = Xp(y)fp{y)- Multiply both sides of 

m = xp{y)fp{y) = ^J2fi-^)xpiy)e'''''^'' 

by e^'^*^^ and integrate over —P/2 < y < P/2. On one hand, 

f{y)e'-^^ydy= / f{y)e'-^^^ydy = f{x), 

p/2 Jr 

and on the other hand. 
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= Enez/(-f)'^(^ + f), 

where 

1 

= / e^"'"^ dy = sinc(7r^). 

P J-P/2 

□ 

4.7 Aliasing 

Obviously most functions are not band limited. When a function is not band 
hmited but the right-hand side of the above "reconstruction formula" is used 

anyway, the error creates an effect called "aliasing." This also occurs when 
one uses the reconstruction formula for a smaple rate lower than the Nyquist 
rate. 

Aliasing is a major concern in the analog-to-digital conversion of video 
and audio signals. 

Theorem 48 (Aliasing Theorem) Assume f is as in the above theorem and 
that the FT of f e L^{R) satisfies 

- (l + li/l)"' 
for some constants ^4 > and a > 0. Then 

fix) = E /(^)smc(2Lx- -n) + E{x), 
where E{x) is an "error" hounded by 
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5 Discrete sine and cosine transforms 



Recall that, given a differentiable, real-valued, periodic function f{x) of pe- 
riod P — 2L, there are a„ with n > and 6„ with n> 1 such that f{x) has 
(real) Fourier series 

FSu{f){x) = Y + 2^[anCOs(-^)+6„sin(-^)]. 

n=l 

where 



and 




When / is even then the 6„ = and we call the FS a cosine series. When 
/ is odd then the a„ = and we call the FS a sine series. In either of these 
cases, it suffices to define / on the interval < x < L instead of < x < P. 

Let us first assume / is even and"discretize" the integral for the k-th 
coefficient of the cosine series and use that as a basis for defining the discrete 
cosine transform or DCT. Using the "left-hand Riemann sum" approximation 
for the integral using N subdivisions, we have 

«fe = l/o^/(^)cos(^)dx 

- I EfJ m/N) cos(^) (#) (23) 

= |Ef=oV(^j7iV)cos(7rA;j7iV). 
This motivates the following definition. 

Definition 49 The N-point discrete cosine transform (or DCT) of the vector 
/ = (/o,...,/iv-i)eM^ is 

N-l 

DCT{f)k=Y,fjCosi7rkj/N), 

3=0 

where <k < N. 
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This transform is represented by the N x N real symmetric matrix 
(cos(7rA;j77V))o<j,fe<Ar-i- 



Example 50 When N — 5, we have 



/ 


1 


1 


1 




1 




1 








1 


1 


1 


1 


1 \ 




1 


cos(f ) 


cos(4i!^) 
cos(Sg-) 


COS 




COS 










1 


0.8090 


0.3090 


-0.3090 


-0.8090 




1 
1 


cos{^) 
cos(Sg-) 


COS 
COS 


f ) 

127r \ 
5 


COS 
COSI 








I 


1 
1 


0.3090 
-0.3090 


-0.8090 
-0.8090 


-0.8090 
0.8089 


0.3090 
0.3090 


V 


1 


cos(^) 


cos(^) 


cos( 


COSI 


5 


J 




1 


-0.8090 


0.3090 


0.3090 


-0.8089 / 



/ 

sm(f) sin(^) sin(^) 

sin(^) Bm(4%^) sin(^) 

sin(^) sin(^) sin(^) 

V sin(^) sin(&) sin(lft) 



\ / 0.0000 0.0000 

Bin(-f^) 0.0000 0.5877 

Bm(^) = 0.0000 0.9610 

Bin(i|^t) 0.0000 0.9610 

sin(i&) / V 0.0000 0.6877 



0.0000 0.0000 0.0000 \ 

0.9610 0.9610 0.6877 

0.6877 -0.6877 -0.9610 

-0.6877 -0.6878 0.9610 

-0.9610 0.9610 -0.6878 / 



Here is the SAGE code for producing the DCT above:: 



sage: RRR = RealField(15) 

sage: MS = MatrixSpace(RRR,5,5) 

sage: r = lambda j: [RRR(cos(pi*j*k/5)) for k in range(5)] 

sage: dct = MS([r(j) for j in range (5)]) 

Next, assume / is odd and "discretize" the integral for the k-th coeffi- 
cient of the sine series and use that as a basis for defining the discrete sine 
transform or DST. Using the "left-hand Riemann sum" approximation for 
the integral using N subdivisions, we have 

h =|/oV(^)sin(^)dx 

^ i EjTo' fiLj/N) sin(^) (^) (24) 

This motivates the following definition. 

Definition 51 The N -point discrete sine transform (or DST) of the vector / = 
(/o = 0,/i,...,/jv-i)gR^ is 

N-l 

DST{f\=J2fjM^kj/N), 

where < k < N. 
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This transform is represented by the N x N real symmetric matrix 
(sin(7rzA;j'/A^))o<j,A:<Af-i- Since the 0-th coordinate is always = 0, since sin(O) = 
0, sometimes this is replaced by a map DST : M^~^ — > represented 
by the {N — 1) x (A^ — 1) real symmetric matrix (sin(7rA;j'/A^))i<j^jfc<Ar-i. 

One difference between these definitions and the definition of the DFT is 
that here the N samples are taken from (0, L), whereas for the DFT the N 
samples are taken from (0, P). 

Is there some way to compute the DOT and the DST in terms of the DFT 
of a function? Let / = (/o, /i, /at, f2N-i) and compute, for < k < N, 

DFT{f), = Y.T=^' f.e-^^^l^ 

= EfJo f^^""'''"" + (-1)' Ef=o' 

= Ef=o'(/. + (-1)'/^+.) cos(7rfcj7iV) 

+^ E7=o'(/. + (-1)'/^+.) sin(7rA;j7iV) 
= DCT{g+ {-lYh)k + iDST{g+ {-lfh)k, 

where g = (/o, /i, /n-i) and h = (/at, f^+i, /27v-i)- Here the DFT is 
based on 2N sample values, whereas the DCT and DST are each based on 
N sample values. 

Let g = {go, gi, ..■,gN-i) G and let g^ = (go, gi, gN-i, 0, 0) G M^^ 
denote its "extension by zero" to M^^. The above computation implies 



DCT{g), = Re{DFT{g,)k), DST{g)u = lm{DFT{g,)k), 

for < A; < A^. In particular, if we can find a "fast" way of computing the 
DFT then we can also compute the DCT and the DST quickly. We turn to 
such efficient computational procedures inthe next section. 



6 Fast Fourier transform 

The fast Fourier transform is a procedure for computing the discrete Fourier 
transform which is especially fast. The term FFT often loosely refers to a 
hybrid combination of the two algorithms presented in this section. 

The algorithm described first, due to James Cooley and John Tukey |CTj . 
works when the number of sample values is a power of 2, say N = 2^, for 
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some integer r > 1. This special case is also referred to as the radix-2 
algorithm. This is the one we will describe in the next section. 



6.1 Cooley-Tukey algorithm (radix-2) 

First, we assume that the powers of W (namely, 1, W , , ^) have 

been precomputed. Note that the computation of the DFT on requires 
iV^ multiplications. This is because the matrix F^r is N x N and each matrix 
entry is involved in the computation of the vector DFT{f) G C^. If M{N) 
denotes the number of multiplications required to compute the DFT then 
the above reasoning shows that 

M{N) < iVl 

To improve on this, the Cooley-Tukey procedure is described next. Let 
N = 2M for the argument below. To be clear about the notation, let = 
^2ni/N^ SO Wl = Wm- Let / = (/o, /i, /iv-i) G be the vector we want 
to compute the DFT of and write 



(25) 



feven — (fo, /2, fN-2) G C*^, fodd — (/l, /s, In-i) G C*^. 

We have, ioi < k < N, 

DFT{f\ =EI-oWn 

- 2^j=o hj vv N + l^j=o hj+i N 

= EfJ hMi +5^ f^^+^^M 

= DFT{l,en)k + w''^DFT{fodd)k. 
Theorem 52 For each N = 2^, M{N) <N -{1 + 2). 

proof: We prove this by mathematical induction on L. This requires 
proving a (base case) step = 4 and a step where we assume the truth of 
the inequality for N/2 and prove it for A^. 

The number of multiplications required to compute the DFT when = 4 
is < 16. Therefore, M(4) < 4 ■ (2 + 2). 

Now assume M{N/2) < f • (L + 1). The Cooley-Tukey procedure (j^ 
shows that M{N) < 2 ■ M{N/2) + A^/2. This and the induction hypothesis 
implies 
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M{N) < 2 ■ M{N/2) + N/2 < 2{— ■ {L + 1)) + iV/2 = (L + 1 + l/2)N. 
This is < A^- (L + 2). □ 

Example 53 When N = 8, the DFT matrix is given by 



( 1 


1 


1 


1 


1 


1 


1 


1 


\ 


1 


Cs 


S8 


S8 


-1 


-Cs 


S8 


-CI 




1 




-1 


S8 


1 


CI 


-1 


-CI 




1 


S8 


S8 


C8 


-1 


-CI 


ss 


-C8 




1 


-1 


1 


-1 


1 


-1 


1 


-1 




1 


-C8 


S8 


-CI 


-1 


Cs 


S8 


S8 




1 


S8 


-1 


CI 


1 


S8 


-1 


S8 






^8 


S8 


-C8 


-1 


ss 


^8 


Cs 


/ 



and the DFT matrix for N/2 



( 1 

1 
1 

Let f = (0,1,2,3,4,5,6,7), so U 
compute 



4 is given by 
1 



1 

C4 
-1 

-C4 



-4C 



28 

■4CI 
-4CI 



1 \ 

-1 -C4 

L -1 

-1 C4 / 

(0,2,4,6), 



(1,3,5,7). We 



4Cs 
4 



4 



-4^ + 4(1-4(8-4 
-4 

4CI-4CI + 4C8-4 

4C| - 4 
V 4CI + 4CI + 4C8-4 / 

k 7„ 



Lei Z denote the diagonal matrix with Q 's on the diagonal, < k < 7 . (This 
matrix represents the factors Wj^ in the formula 112^) above.) It turns out to 
be more useful for computations to split this matrix up into two parts: 
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z. = 



Note Zo 



/ 1 \ 

C8 
CI 
Vo CI/ 

'Zi. For the first 4 coordinates, 



( 



F^feven + ^'1^4/, 



odd 





/CI 






















V 





we use 






28 






4CI- 


4C8 




4CI- 


4 


+ 4CI - 


4C8 



and for the last 4 coordinates, we use 

/ -4 \ 

Ff ,7Ff - 4CI-4CI + 4C8-4 

^ijeven "T ^2^Ajodd — ^^2 _ ^ 

V4CI + 4CI + 4C8-4/ 
We see that, as (H^j predicts, Fgf is equal to [Fifeven + ZiF^fodd, F^feven + 

Z2F4:fodd\- 

Here is the SAGE code for the above computations. 



sage : 


MS4 


= MatrixSpace(CyclotomicField(4) ,4,4) 




sage : 


MS8 


= MatrixSpace(CyclotomicField(8) ,8,8) 




sage : 


V4 = 


VectorSpace (CyclotomicField(4) ,4) 




sage : 


V8 = 


VectorSpace (CyclotomicField(8) ,8) 




sage : 


z4 = 


CyclotomicField(4) .genO 




sage : 


z8 = 


CyclotomicField(8) .genO 




sage : 


r4 = 


lambda k: [z4~(j*k) for j in range (4)] 




sage : 


r8 = 


lambda k: [z8~(j*k) for j in range (8)] 




sage : 


F4 = 


MS4([r4(k) for k in range (4)]) 




sage : 


F8 = 


MS8([r8(k) for k in range (8)]) 




sage : 


f = V8([0,l,2,3,4,5,6,7]) 




sage : 


fe = 


V4([0,2,4,6]) 




sage : 


fo = 


V4([l,3,5,7]) 




sage : 


FFTe 


= [(F4*fe) [j]+z8"j*(F4*fo) [j] for j in 


range (4)] 


sage : 


FFTo 


= [(F4*fe) [j]-z8"j*(F4*fo) [j] for j in 


range (4)] 


sage : 


FFTe+FFTo 




[28, 









-4*zeta8~3 
-4*zeta8~2 



4*zeta8"2 
4, 



4*zeta8 - 4, 
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-4*zeta8~3 + 4*zeta8~2 - 4*zeta8 - 4, 
-4, 

4*zeta8~3 - 4*zeta8'2 + 4*zeta8 - 4, 
4*zeta8~2 - 4, 

4*zeta8~3 + 4*zeta8'2 + 4*zeta8 - 4] 
sage: [(F8*f ) [j] for j in range(8)] 

[28, 

-4*zeta8~3 - 4*zeta8"2 - 4*zeta8 - 4, 
-4*zeta8~2 - 4, 

-4*zeta8~3 + 4*zeta8"2 - 4*zeta8 - 4, 
-4, 

4*zeta8~3 - 4*zeta8'2 + 4*zeta8 - 4, 
4*zeta8~2 - 4, 

4*zeta8~3 + 4*zeta8'2 + 4*zeta8 - 4] 



Finally, we give an example which only illustrates SAGE 's implementation 
of the FFT (which calls functions in the GSL |GSLj ). as compared to its 
implementation of its DFT (which is implemented in Python but calls Pari 
for the computations involving A^-th roots of unity over Q): 

sage: J = range (30) 

sage: A = [QQ (int (10* (randomO -1/2) ) ) for i in J] 
sage: s = IndexedSequence(A, J) 
sage: time dfts = s.dftO 

CPU times: user 0.86 s, sys: 0.04 s, total: 0.90 s 

Wall time: 0.94 

sage: time ffts = s.fftO 

CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s 
Wall time: 0.00 
sage: J = range (3000) 

sage: A = [QQ (int (10* (randomO -1/2) ) ) for i in J] 
sage: s = lndexedSequence(A, J) 
sage: time ffts = s.fftO 

CPU times: user 0.21 s, sys: 0.00 s, total: 0.21 s 
Wall time: 0.21 

As you can see, for a sample vector in C^*"^", SAGE can compute the FFT in 
about |-th of a second. However, of you try to compute the DFT using this 
example, SAGE will probably give you an error related to its extreme size. 
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6.2 Rader's algorithm 

In this subsection, we assume N is prime. Here we breifly describe an algo- 
rithm due to Rader |Raj for computing the DFT on Vn for prime N. 

The basic idea is to rewrite the DFT on as a convolution on Vat-i- 
Remark n is then used to show that this convolution can be computed using 
a "fast" algorithm. 

The first step in the algorithm is to select an element g, 1 < g < N — 1, 
which has the property that every element y in {1,2, N — 1} can be written 
in the form y = g^ (mod A^) for some x. This element g s called a primitive 
root (mod N) (or a generator of {Z/NZ^), where (Z/A^Z)^ = Z/NZ - 



{0}. Here is 


a 


table of primitive roots for 


various small primes A^, and 


a 


demonstration, 


in the 


case = 17, that g = 


= 3 is indeed a primitive root: 








P 


3 5 7 11 11 13 


17 19 23 








9 


2 2 3 2 2 2 


3 2 5 




k 





1 2 


3 4 5 6 7 


8 9 10 11 12 13 14 


15 


g" (mod N) 


1 


3 9 


10 13 5 15 11 


16 14 8 7 4 12 2 


6 



The SAGE command which produces the smallest primitive root (mod A^) 
is primitive_root(N). For example, the above tables were produced using 
the commands 



sage: [priniitive_root (p) for p in [3,5,7,11,13,17,19,23]] 

^2 J 2j 3j 2j 2j 3j 2j 5^ 
sage: N = 17; g = 3 
sage: [g~k°/oN for k in range (N-1)] 

[1, 3, 9, 10, 13, 5, 15, 11, 16, 14, 8, 7, 4, 12, 2, 6] 

The next step in the algorithm is to rewrite the DFT on Vn as a convo- 
lution on Vn-i- To this end, fix / e Vn and define hi, /12 € Vat-i as 

h{x) = f{g-n, /i2(x)=e-^-^^/^, 
for X G {0, 1, 2, A^ - 2} = Z/(A^ - 1)Z. For = we compute 

DFT{f)k = DFT{f)o = /(O) + /(I) + ... + f{N - 1). 
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For £ 7^ 0, we let m = m{£) denote the element of {0, 1, 2, A^ — 2} such that 
£ = g~"^. For 7^ 0, we let n = n{k) denote the element of {0, 1,2, A^ — 2} 
such that k = g^. Now write 

DFT{f), = Ef=oV(^)e-^^^'=^/^ 

= /(O) + Em=o hi{m)h2{n - m) 
= f{Q) + h^*h2{n). 

This is a convolution on Vn-i- If iV — 1 is a power of 2 (e.g., for = 17) 
then use Remark ^ and the radix-2 Cooley-Tukey algorithm described in the 
previous section to quickly compute hi * If — 1 is not a power of 2, 
the best thing to do is to let P denote the smallest power of 2 greater than 
A^ — 1 and extend both hi and /i2 by to the range {0, 1, 2, P — 1} (this is 
called "padding" the functions). Call these new functions hi and /i2. Using 
DFT{f)k = /(O) + hi * h2{n) (where k = g'"-). We now have expressed the 
DFT on Vm in terms of a convolution on Vp, where P < 2N. By Remark 
^ we know that this can be computed in < cA^log2(A^) multiplications, for 
some constant c > 1 (in fact, c = 4 should work if A^ is sufficiently large). 

7 Fourier optics 

Fourier optics is a special topic in the theory of optics. Good references are 
Goodman \Go\ and chapter 5 of Walker jWlj . This section owes much to 
discussions with Prof. Larry Tankersley of the USNA Physics dept. 

The object of this section is to describe a method for computing cer- 
tain quantities arising in defraction experiments. I'll try to describe these 
experiments next. 

Consider a monochromatic light source having wavelength A. For ex- 
ample, the light emitted from a laser would fit thi^ Imagine x- and y- 
coordinates in the aperture plane. The aperture function A{x,y) is defined 

^In reality, a laser beam is too narrow, so a series of lenses is required to widen it and 
then straighten out the light rays. You need to make sure that the beam is wide enough 
for it to be possible to place a small aperture (a slit or defraction grating, for instance) in 
front of it to allow only the light through that the aperture. 
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to be 1 if light can pass though the "sht" at (x, y) and otherwisJ^ The 
hght which passes through the aperture is pictured on a screen. It is this 
pattern which we wish to describe in this section. 



ight rays 




aperature plane 



Figure 29: Sht experiment set-up for Fourier optics modeL 



^°(It is even possible to image values of A{x, y) in the range between and 1 representing 
a partially opaque screen, though wc shall not need this here. However, we do assume A 
is real- valued. 
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When the aperture is a square (whose sides are aUgned parallel to the x- 
and y-axes), the image projected on the screen of this experiment looks like 
a dotted plus sign, where the dots get fainter and fainter as they move away 
from the center. A square slit diffraction pattern is pictured belo\\0 







• 






It 







Figure 30: Square aperture diffraction experiment, Betzler |B] . 

The goal of this last section will be to describe the mathematics behind this 
"dotted plus sign" square slit diffraction pattern. 

7.1 The mathematical model 

The theory we shall sketch is called scalar defraction theory. A special case 
which we shall concentrate on is the Fraunhofer defraction model. 

Let L denote the distance from the aperture to the screen, A (as above) 
the wavelength of the light, and a > the width of the slit. The Fresnel 
numbed is the quantity F = When F > 1 the screen is "sufficiently 
far" (past the "Fresnel threshold") from the aperture that the wavefronts 
created by the slit have negligable curvature. 



^^This image is copyright Prof. Dr. Klaus Bctzlcr and reproduced with his kind per- 
mission (stated in an email dated May 2, 2007) 
^^Pronounced "Fre-nell". 
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Here is Wikipedia's imag^ |WFdj : 



4 




u ■ u 



p 



Figure 31: On this diagram, a wave is diffracted and observed at point a. As this 
point is moved further back, beyond the Fresnel threshold, Fraunhofer diffraction 
occurs. 



We also assume that the index of refraction of the medium is 1. If a light 
ray of wavelength A travels from P to Q, points at a distance of r = | |Q — P| | 
from each other, at time t then the light amplitude ijj{Q,t) at Q arising from 
this light ray satisfies the property 



where ip{P,t) is the light amplitude at P. We can and do assume P = 
(x, y, 0) is a point in the aperture plane and Q = {x', y', L) is a point on the 
screen. We assume ||Q — P|| is not too big - so that light travels essentially 
instantaneously from P to Q. The superposition principle implies that the 
total light amplitude at Q satisfies 



where we have identified the plane of the aperture slit with the Cartesian 
plane M?. If the distance between the aperture and he screen is large enough 
then r is a constant 

The light intensity is defined by 

^•^Wikipedia's images are copyright the Wikipedia Foundation and distributed under 
the GNU Documentation License. 



^(Q,t)=7A(P,t) 



Ar ' 




27ri||Q-P|| 



(26) 
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/(g)= lim i rmQ,t)\'dt. 

The light amphtude is not typically measurable (at least not in air, with 
present day technology), but the intensity is. 

Let P and P' be points on the aperture plane where A{P) ^ and 
A{P') 7^ 0. The coherency function is defined by 

T(P,P') = lim , ^ ^ [ ip(P,tU(P',t) dt. 
^ ' ' ^I{P)I{P') T Jo ' ' ' 

We say thet the light is coherent if T{P, P') = 1 for all such P, P' . 
Using (f^ . we have 

J(Q) = limT^oo ^ ij{Q,t)ij{Q,t) dt 

= hm^^^ i A{P)i,{P, t)^=p^ dP) X 

X (/j,, A{P')W^Y^^^ dP') dt 

_ . r A(P)A{P') 2..(||Q-P|| IIQ-P'II) 

— mny^oo Jk2 J]r2 a2||Q -P||||Q-P '||^ ^ 

X /o^ V'(^, t)ij{P',t)dP dP') dt. 
Assuming that the light is coherent, this is 



_ r r A(P)A{P') 2..(||q-p|| ||q-p'||) ^-^^-^ 

We assumJ^ that /(-P) = 1 for all points P on the aperture screen with 
A{P) 7^ 0. In this case, the above reasoning leads to 

This is the key equation in scalar diffraction theory that enables one to 
approximate the intensity function in terms of the Fourier transform of the 
aperture function. 

^"^The assumption that the intensity is constant at the aperture is all that is really 
necessary. We assume that this constant is = 1 for simplicity. 
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7.2 The Fraunhofer model 



Notation: In x, 7/, ^-coordinates, the aperature plane is described hy z = 0, 
the screen is described hj z = L > 0. In the diagram below, if Q — {x, y, L) 
then Q' = {x', y\ L), P' = {x, y, 0), and P = {x', y\ 0). 



Q' 
Q 



aperature plane 



Figure 32: Notation for slit experiment set-up for Fraunhofer model. 

In particular, the vector v — Q—P' is orthogonal to the vector w — P' — P 
and P ■ P' — P ■ Q. These give us 



\Q-P\\^\\Q-P' + p'-P\\ = (||Q_p'||2 + ||p/_p||2)l/2 



\P'-P\\^\l/2 T , IIP'-PII 



1,2 

1 , 



2L 



+ ... , 



by the power series expansion (1 + x)^/^ = 1 + |x + .... In addition to the 
coherence assumption, we also assume that L is so large and the aperture 
opening (i.e., the support of the aperture function A in the aperture plane, 
which we identify with M^) is so small that the error in 



27ri||Q-P| 7ri||P'-P|p 2TTiL 



'e ^ 
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is negligable. We expand \ \P' - P\Y = | |P'| p - 2P' ■ P + | |P| p, so 

e A e Lx e \ e L\ e i-^' = e ^-^ e ^ e e . 

If / G define the 2- dimensional Fourier transforms by 



We also assume that L is so large and the aperture opening is so small that 
the error in e L>r~ ^ 1 is negligable. Plugging this into (jTTj) . therefore gives 

In otherwords, if Q = {x,y,L) then /(Q) = px^ where A de- 

notes the 2-d Fourier transform. 

Example 54 VFe consider the example where the slit is a small rectangle. 
For example, suppose 

A{x,y) = X-ei,e^{x)X-e2,e2{y), 

where ei > and €2 > are given, and Xa,b represents the function which is 
1 for a < X << b and 0, otherwise. 
We know 



therefore 



and so 



J-a ™W 



m) = ^(sinc(2vrei-^)sinc(2vre2-^))2. 
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Here are the SAGE commands which produce the plot below for this in- 
tensity function (the axes have been scaled for simplicity). 

sage: f = " (sin(x) /x) '2* (sin(y) /y) '2" 
sage: opts = " [plot_f ormat , openmath] " 

sage: maxima. plotSd (f, "[x, -5, 5]", "[y, -5, 5]", opts) 




Figure 33: Intensity function for a square slit experiment in the Fraunhofer 
model. 

This plot is consistent with the image pictured in the '^dotted plus sign" square 
slit diffraction pattern. As the mathematical explanation of this image was 
the goal of this last section of these notes, we are done. 

8 Additional reading 

For an elementary introduction to wavelets, appropriate for capstone projects, 
see Frazier jF]. 
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A number of papers from the American Math. Monthly in tlie theory of 
FS is available at the URL 

*W;tp : //math, fullerton . edu/mathews/ c2Q03/FourierTransf ormBib/Liiiks/FourierTransf o^B^^^k^2Tht^^ 

A number of those articles would also make good starters for a capstone 
project. 

Finally, we recommend Walker |W2j and Korner jKj as excellent addi- 
tional references. 
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